Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 971:e54c2f54dbfe feature/getBoundaryOperator

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic. Use only the changes made to multiblock.DiffOp and scheme.Elastic2dVariable. DiffOp.getBoundaryOperator/Quadrature now use scheme methods instead of propeties.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 25 Dec 2018 07:50:07 +0100
parents 23d9ca6755be
children 104f0af001e0
comparison
equal deleted inserted replaced
957:2412f407749a 971:e54c2f54dbfe
28 28
29 % Traction operators used for BC 29 % Traction operators used for BC
30 T_l, T_r 30 T_l, T_r
31 tau_l, tau_r 31 tau_l, tau_r
32 32
33 H, Hi % Inner products 33 H, Hi, H_1D % Inner products
34
35 phi % Borrowing constant for (d1 - e^T*D1) from R
36 gamma % Borrowing constant for d1 from M
37 H11 % First element of H
38
39 % Borrowing from H, M, and R
40 thH
41 thM
42 thR
43
44 e_l, e_r 34 e_l, e_r
35 e_w, e_e, e_s, e_n
36
37
45 d1_l, d1_r % Normal derivatives at the boundary 38 d1_l, d1_r % Normal derivatives at the boundary
46 E % E{i}^T picks out component i 39 E % E{i}^T picks out component i
47 40
48 H_boundary % Boundary inner products 41 H_boundary % Boundary inner products
42 H_w, H_e, H_s, H_n
49 43
50 % Kroneckered norms and coefficients 44 % Kroneckered norms and coefficients
51 RHOi_kron 45 RHOi_kron
52 Hi_kron 46 Hi_kron
47
48 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
49 theta_R % Borrowing (d1- D1)^2 from R
50 theta_H % First entry in norm matrix
51 theta_M % Borrowing d1^2 from M.
52
53 % Structures used for adjoint optimization
54 B
53 end 55 end
54 56
55 methods 57 methods
56 58
57 function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) 59 % The coefficients can either be function handles or grid functions
60 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
58 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); 61 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
59 default_arg('lambda_fun', @(x,y) 0*x+1); 62 default_arg('lambda', @(x,y) 0*x+1);
60 default_arg('mu_fun', @(x,y) 0*x+1); 63 default_arg('mu', @(x,y) 0*x+1);
61 default_arg('rho_fun', @(x,y) 0*x+1); 64 default_arg('rho', @(x,y) 0*x+1);
62 dim = 2; 65 dim = 2;
63 66
64 assert(isa(g, 'grid.Cartesian')) 67 assert(isa(g, 'grid.Cartesian'))
65 68
66 lambda = grid.evalOn(g, lambda_fun); 69 if isa(lambda, 'function_handle')
67 mu = grid.evalOn(g, mu_fun); 70 lambda = grid.evalOn(g, lambda);
68 rho = grid.evalOn(g, rho_fun); 71 end
72 if isa(mu, 'function_handle')
73 mu = grid.evalOn(g, mu);
74 end
75 if isa(rho, 'function_handle')
76 rho = grid.evalOn(g, rho);
77 end
78
69 m = g.size(); 79 m = g.size();
70 m_tot = g.N(); 80 m_tot = g.N();
71 81
72 h = g.scaling(); 82 h = g.scaling();
73 lim = g.lim; 83 lim = g.lim;
85 ops{i} = opSet{i}(m(i), lim{i}, order); 95 ops{i} = opSet{i}(m(i), lim{i}, order);
86 end 96 end
87 97
88 % Borrowing constants 98 % Borrowing constants
89 for i = 1:dim 99 for i = 1:dim
90 beta = ops{i}.borrowing.R.delta_D; 100 obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D;
91 obj.H11{i} = ops{i}.borrowing.H11; 101 obj.theta_H{i} = h(i)*ops{i}.borrowing.H11;
92 obj.phi{i} = beta/obj.H11{i}; 102 obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1;
93 obj.gamma{i} = ops{i}.borrowing.M.d1;
94
95 % Better names
96 obj.thR{i} = ops{i}.borrowing.R.delta_D;
97 obj.thM{i} = ops{i}.borrowing.M.d1;
98 obj.thH{i} = ops{i}.borrowing.H11;
99 end 103 end
100 104
101 I = cell(dim,1); 105 I = cell(dim,1);
102 D1 = cell(dim,1); 106 D1 = cell(dim,1);
103 D2 = cell(dim,1); 107 D2 = cell(dim,1);
144 % Boundary operators 148 % Boundary operators
145 obj.e_l{1} = kron(e_l{1},I{2}); 149 obj.e_l{1} = kron(e_l{1},I{2});
146 obj.e_l{2} = kron(I{1},e_l{2}); 150 obj.e_l{2} = kron(I{1},e_l{2});
147 obj.e_r{1} = kron(e_r{1},I{2}); 151 obj.e_r{1} = kron(e_r{1},I{2});
148 obj.e_r{2} = kron(I{1},e_r{2}); 152 obj.e_r{2} = kron(I{1},e_r{2});
153 obj.e_w = obj.e_l{1};
154 obj.e_e = obj.e_r{1};
155 obj.e_s = obj.e_l{2};
156 obj.e_n = obj.e_r{2};
149 157
150 obj.d1_l{1} = kron(d1_l{1},I{2}); 158 obj.d1_l{1} = kron(d1_l{1},I{2});
151 obj.d1_l{2} = kron(I{1},d1_l{2}); 159 obj.d1_l{2} = kron(I{1},d1_l{2});
152 obj.d1_r{1} = kron(d1_r{1},I{2}); 160 obj.d1_r{1} = kron(d1_r{1},I{2});
153 obj.d1_r{2} = kron(I{1},d1_r{2}); 161 obj.d1_r{2} = kron(I{1},d1_r{2});
181 obj.H = kron(H{1},H{2}); 189 obj.H = kron(H{1},H{2});
182 obj.Hi = inv(obj.H); 190 obj.Hi = inv(obj.H);
183 obj.H_boundary = cell(dim,1); 191 obj.H_boundary = cell(dim,1);
184 obj.H_boundary{1} = H{2}; 192 obj.H_boundary{1} = H{2};
185 obj.H_boundary{2} = H{1}; 193 obj.H_boundary{2} = H{1};
194 obj.H_1D = {H{1}, H{2}};
195 obj.H_w = H{2};
196 obj.H_e = H{2};
197 obj.H_s = H{1};
198 obj.H_n = H{1};
186 199
187 % E{i}^T picks out component i. 200 % E{i}^T picks out component i.
188 E = cell(dim,1); 201 E = cell(dim,1);
189 I = speye(m_tot,m_tot); 202 I = speye(m_tot,m_tot);
190 for i = 1:dim 203 for i = 1:dim
211 D2_mu{j}*E{i}' ... 224 D2_mu{j}*E{i}' ...
212 ); 225 );
213 end 226 end
214 end 227 end
215 obj.D = D; 228 obj.D = D;
216 %=========================================% 229 %=========================================%'
217 230
218 % Numerical traction operators for BC. 231 % Numerical traction operators for BC.
219 % Because d1 =/= e0^T*D1, the numerical tractions are different 232 % Because d1 =/= e0^T*D1, the numerical tractions are different
220 % at every boundary. 233 % at every boundary.
221 T_l = cell(dim,1); 234 T_l = cell(dim,1);
235 T_l{j} = cell(dim,dim); 248 T_l{j} = cell(dim,dim);
236 T_r{j} = cell(dim,dim); 249 T_r{j} = cell(dim,dim);
237 tau_l{j} = cell(dim,1); 250 tau_l{j} = cell(dim,1);
238 tau_r{j} = cell(dim,1); 251 tau_r{j} = cell(dim,1);
239 252
253 LAMBDA_l = e_l{j}'*LAMBDA*e_l{j};
254 LAMBDA_r = e_r{j}'*LAMBDA*e_r{j};
255 MU_l = e_l{j}'*MU*e_l{j};
256 MU_r = e_r{j}'*MU*e_r{j};
257
258 [~, n_l] = size(e_l{j});
259 [~, n_r] = size(e_r{j});
260
240 % Loop over components 261 % Loop over components
241 for i = 1:dim 262 for i = 1:dim
242 tau_l{j}{i} = sparse(m_tot,dim*m_tot); 263 tau_l{j}{i} = sparse(n_l, dim*m_tot);
243 tau_r{j}{i} = sparse(m_tot,dim*m_tot); 264 tau_r{j}{i} = sparse(n_r, dim*m_tot);
244 for k = 1:dim 265 for k = 1:dim
245 T_l{j}{i,k} = ... 266 T_l{j}{i,k} = ...
246 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... 267 -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})...
247 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... 268 -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})...
248 -d(i,k)*MU*e_l{j}*d1_l{j}'; 269 -d(i,k)*MU_l*d1_l{j}';
249 270
250 T_r{j}{i,k} = ... 271 T_r{j}{i,k} = ...
251 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... 272 d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})...
252 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... 273 +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})...
253 +d(i,k)*MU*e_r{j}*d1_r{j}'; 274 +d(i,k)*MU_r*d1_r{j}';
254 275
255 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; 276 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
256 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; 277 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
257 end 278 end
258 279
259 end 280 end
260 end 281 end
282
283 % Transpose T and tau to match boundary operator convention
284 for i = 1:dim
285 for j = 1:dim
286 tau_l{i}{j} = transpose(tau_l{i}{j});
287 tau_r{i}{j} = transpose(tau_r{i}{j});
288 for k = 1:dim
289 T_l{i}{j,k} = transpose(T_l{i}{j,k});
290 T_r{i}{j,k} = transpose(T_r{i}{j,k});
291 end
292 end
293 end
294
261 obj.T_l = T_l; 295 obj.T_l = T_l;
262 obj.T_r = T_r; 296 obj.T_r = T_r;
263 obj.tau_l = tau_l; 297 obj.tau_l = tau_l;
264 obj.tau_r = tau_r; 298 obj.tau_r = tau_r;
265 299
272 obj.m = m; 306 obj.m = m;
273 obj.h = h; 307 obj.h = h;
274 obj.order = order; 308 obj.order = order;
275 obj.grid = g; 309 obj.grid = g;
276 obj.dim = dim; 310 obj.dim = dim;
311
312 % B, used for adjoint optimization
313 B = cell(dim, 1);
314 for i = 1:dim
315 B{i} = cell(m_tot, 1);
316 end
317
318 for i = 1:dim
319 for j = 1:m_tot
320 B{i}{j} = sparse(m_tot, m_tot);
321 end
322 end
323
324 ind = grid.funcToMatrix(g, 1:m_tot);
325
326 % Direction 1
327 for k = 1:m(1)
328 c = sparse(m(1),1);
329 c(k) = 1;
330 [~, B_1D] = ops{1}.D2(c);
331 for l = 1:m(2)
332 p = ind(:,l);
333 B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
334 end
335 end
336
337 % Direction 2
338 for k = 1:m(2)
339 c = sparse(m(2),1);
340 c(k) = 1;
341 [~, B_1D] = ops{2}.D2(c);
342 for l = 1:m(1)
343 p = ind(l,:);
344 B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
345 end
346 end
347
348 obj.B = B;
277 349
278 end 350 end
279 351
280 352
281 % Closure functions return the operators applied to the own domain to close the boundary 353 % Closure functions return the operators applied to the own domain to close the boundary
293 comp = bc{1}; 365 comp = bc{1};
294 type = bc{2}; 366 type = bc{2};
295 367
296 % j is the coordinate direction of the boundary 368 % j is the coordinate direction of the boundary
297 j = obj.get_boundary_number(boundary); 369 j = obj.get_boundary_number(boundary);
298 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 370 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
371
299 372
300 E = obj.E; 373 E = obj.E;
301 Hi = obj.Hi; 374 Hi = obj.Hi;
302 LAMBDA = obj.LAMBDA; 375 LAMBDA = obj.LAMBDA;
303 MU = obj.MU; 376 MU = obj.MU;
314 switch type 387 switch type
315 388
316 % Dirichlet boundary condition 389 % Dirichlet boundary condition
317 case {'D','d','dirichlet','Dirichlet'} 390 case {'D','d','dirichlet','Dirichlet'}
318 391
319 phi = obj.phi{j}; 392 alpha = obj.getBoundaryOperator('alpha', boundary);
320 h = obj.h(j);
321 h11 = obj.H11{j}*h;
322 gamma = obj.gamma{j};
323
324 a_lambda = dim/h11 + 1/(h11*phi);
325 a_mu_i = 2/(gamma*h);
326 a_mu_ij = 2/h11 + 1/(h11*phi);
327
328 d = @kroneckerDelta; % Kronecker delta
329 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
330 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
331 + d(i,j)* a_mu_i*MU ...
332 + db(i,j)*a_mu_ij*MU );
333 393
334 % Loop over components that Dirichlet penalties end up on 394 % Loop over components that Dirichlet penalties end up on
335 for i = 1:dim 395 for i = 1:dim
336 C = T{k,i}; 396 C = transpose(T{k,i});
337 A = -d(i,k)*alpha(i,j); 397 A = -e*transpose(alpha{i,k});
338 B = A + C; 398 B = A + e*C;
339 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); 399 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
340 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 400 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
341 end 401 end
342 402
343 % Free boundary condition 403 % Free boundary condition
344 case {'F','f','Free','free','traction','Traction','t','T'} 404 case {'F','f','Free','free','traction','Traction','t','T'}
345 closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); 405 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}';
346 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; 406 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
347 407
348 % Unknown boundary condition 408 % Unknown boundary condition
349 otherwise 409 otherwise
350 error('No such boundary condition: type = %s',type); 410 error('No such boundary condition: type = %s',type);
351 end 411 end
352 end 412 end
353 413
354 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 414 % type Struct that specifies the interface coupling.
415 % Fields:
416 % -- tuning: penalty strength, defaults to 1.2
417 % -- interpolation: type of interpolation, default 'none'
418 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
419
420 defaultType.tuning = 1.2;
421 defaultType.interpolation = 'none';
422 default_struct('type', defaultType);
423
424 switch type.interpolation
425 case {'none', ''}
426 [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type);
427 case {'op','OP'}
428 [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type);
429 otherwise
430 error('Unknown type of interpolation: %s ', type.interpolation);
431 end
432 end
433
434 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
435 tuning = type.tuning;
436
355 % u denotes the solution in the own domain 437 % u denotes the solution in the own domain
356 % v denotes the solution in the neighbour domain 438 % v denotes the solution in the neighbour domain
357 % Operators without subscripts are from the own domain. 439 % Operators without subscripts are from the own domain.
358 tuning = 1.2;
359 440
360 % j is the coordinate direction of the boundary 441 % j is the coordinate direction of the boundary
361 j = obj.get_boundary_number(boundary); 442 j = obj.get_boundary_number(boundary);
362 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary); 443 j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
363 444
364 % Get boundary operators 445 % Get boundary operators
365 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 446 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
366 [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary); 447 [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_boundary);
367 448
368 % Operators and quantities that correspond to the own domain only 449 % Operators and quantities that correspond to the own domain only
369 Hi = obj.Hi; 450 Hi = obj.Hi;
370 RHOi = obj.RHOi; 451 RHOi = obj.RHOi;
371 dim = obj.dim; 452 dim = obj.dim;
385 lambda_v = e_v'*LAMBDA_v*e_v; 466 lambda_v = e_v'*LAMBDA_v*e_v;
386 mu_v = e_v'*MU_v*e_v; 467 mu_v = e_v'*MU_v*e_v;
387 %------------------------- 468 %-------------------------
388 469
389 % Borrowing constants 470 % Borrowing constants
390 h_u = obj.h(j); 471 theta_R_u = obj.theta_R{j};
391 thR_u = obj.thR{j}*h_u; 472 theta_H_u = obj.theta_H{j};
392 thM_u = obj.thM{j}*h_u; 473 theta_M_u = obj.theta_M{j};
393 thH_u = obj.thH{j}*h_u; 474
394 475 theta_R_v = neighbour_scheme.theta_R{j_v};
395 h_v = neighbour_scheme.h(j_v); 476 theta_H_v = neighbour_scheme.theta_H{j_v};
396 thR_v = neighbour_scheme.thR{j_v}*h_v; 477 theta_M_v = neighbour_scheme.theta_M{j_v};
397 thH_v = neighbour_scheme.thH{j_v}*h_v; 478
398 thM_v = neighbour_scheme.thM{j_v}*h_v; 479 function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu)
399 480 alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M);
400 % alpha = penalty strength for normal component, beta for tangential 481 alpha_ij = mu/(2*th_H) + mu/(4*th_R);
401 alpha_u = dim*lambda_u/(4*thH_u) + lambda_u/(4*thR_u) + mu_u/(2*thM_u); 482 end
402 alpha_v = dim*lambda_v/(4*thH_v) + lambda_v/(4*thR_v) + mu_v/(2*thM_v); 483
403 beta_u = mu_u/(2*thH_u) + mu_u/(4*thR_u); 484 [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u);
404 beta_v = mu_v/(2*thH_v) + mu_v/(4*thR_v); 485 [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v);
405 alpha = alpha_u + alpha_v; 486 sigma_ii = alpha_ii_u + alpha_ii_v;
406 beta = beta_u + beta_v; 487 sigma_ij = alpha_ij_u + alpha_ij_v;
407 488
408 d = @kroneckerDelta; % Kronecker delta 489 d = @kroneckerDelta; % Kronecker delta
409 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 490 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
410 strength = @(i,j) tuning*(d(i,j)*alpha + db(i,j)*beta); 491 sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
411 492
412 % Preallocate 493 % Preallocate
413 closure = sparse(dim*m_tot_u, dim*m_tot_u); 494 closure = sparse(dim*m_tot_u, dim*m_tot_u);
414 penalty = sparse(dim*m_tot_u, dim*m_tot_v); 495 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
415 496
416 % Loop over components that penalties end up on 497 % Loop over components that penalties end up on
417 for i = 1:dim 498 for i = 1:dim
418 closure = closure - E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e'*E{i}'; 499 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
419 penalty = penalty + E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e_v'*E_v{i}'; 500 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
420 501
421 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}; 502 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau{i}';
422 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}; 503 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau_v{i}';
423 504
424 % Loop over components that we have interface conditions on 505 % Loop over components that we have interface conditions on
425 for k = 1:dim 506 for k = 1:dim
426 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}'; 507 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}';
427 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}'; 508 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}';
428 end 509 end
429 end 510 end
511 end
512
513 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
514 error('Non-conforming interfaces not implemented yet.');
430 end 515 end
431 516
432 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 517 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
433 function [j, nj] = get_boundary_number(obj, boundary) 518 function [j, nj] = get_boundary_number(obj, boundary)
434 519
449 end 534 end
450 end 535 end
451 536
452 % Returns the boundary operator op for the boundary specified by the string boundary. 537 % Returns the boundary operator op for the boundary specified by the string boundary.
453 % op: may be a cell array of strings 538 % op: may be a cell array of strings
454 function [varargout] = get_boundary_operator(obj, op, boundary) 539 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
540 function [varargout] = getBoundaryOperator(obj, op, boundary)
455 541
456 switch boundary 542 switch boundary
457 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 543 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
458 j = 1; 544 j = 1;
459 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 545 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
464 550
465 if ~iscell(op) 551 if ~iscell(op)
466 op = {op}; 552 op = {op};
467 end 553 end
468 554
469 for i = 1:length(op) 555 for k = 1:length(op)
470 switch op{i} 556 switch op{k}
471 case 'e' 557 case 'e'
472 switch boundary 558 switch boundary
473 case {'w','W','west','West','s','S','south','South'} 559 case {'w','W','west','West','s','S','south','South'}
474 varargout{i} = obj.e_l{j}; 560 varargout{k} = obj.e_l{j};
475 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 561 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
476 varargout{i} = obj.e_r{j}; 562 varargout{k} = obj.e_r{j};
477 end 563 end
564
565 case 'e_tot'
566 e = obj.getBoundaryOperator('e', boundary);
567 I_dim = speye(obj.dim, obj.dim);
568 varargout{k} = kron(e, I_dim);
569
478 case 'd' 570 case 'd'
479 switch boundary 571 switch boundary
480 case {'w','W','west','West','s','S','south','South'} 572 case {'w','W','west','West','s','S','south','South'}
481 varargout{i} = obj.d1_l{j}; 573 varargout{k} = obj.d1_l{j};
482 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 574 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
483 varargout{i} = obj.d1_r{j}; 575 varargout{k} = obj.d1_r{j};
484 end 576 end
485 case 'H' 577
486 varargout{i} = obj.H_boundary{j};
487 case 'T' 578 case 'T'
488 switch boundary 579 switch boundary
489 case {'w','W','west','West','s','S','south','South'} 580 case {'w','W','west','West','s','S','south','South'}
490 varargout{i} = obj.T_l{j}; 581 varargout{k} = obj.T_l{j};
491 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 582 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
492 varargout{i} = obj.T_r{j}; 583 varargout{k} = obj.T_r{j};
493 end 584 end
585
494 case 'tau' 586 case 'tau'
495 switch boundary 587 switch boundary
496 case {'w','W','west','West','s','S','south','South'} 588 case {'w','W','west','West','s','S','south','South'}
497 varargout{i} = obj.tau_l{j}; 589 varargout{k} = obj.tau_l{j};
498 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 590 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
499 varargout{i} = obj.tau_r{j}; 591 varargout{k} = obj.tau_r{j};
500 end 592 end
593
594 case 'tau_tot'
595 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
596
597 I_dim = speye(obj.dim, obj.dim);
598 e_tot = kron(e, I_dim);
599 E = obj.E;
600 tau_tot = (e_tot'*E{1}*e*tau{1}')';
601 for i = 2:obj.dim
602 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
603 end
604 varargout{k} = tau_tot;
605
606 case 'H'
607 varargout{k} = obj.H_boundary{j};
608
609 case 'alpha'
610 % alpha = alpha(i,j) is the penalty strength for displacement BC.
611 e = obj.getBoundaryOperator('e', boundary);
612
613 tuning = 1.2;
614 LAMBDA = obj.LAMBDA;
615 MU = obj.MU;
616
617 dim = obj.dim;
618 theta_R = obj.theta_R{j};
619 theta_H = obj.theta_H{j};
620 theta_M = obj.theta_M{j};
621
622 a_lambda = dim/theta_H + 1/theta_R;
623 a_mu_i = 2/theta_M;
624 a_mu_ij = 2/theta_H + 1/theta_R;
625
626 d = @kroneckerDelta; % Kronecker delta
627 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
628 alpha = cell(obj.dim, obj.dim);
629
630 alpha_func = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
631 + d(i,j)* a_mu_i*MU ...
632 + db(i,j)*a_mu_ij*MU );
633 for i = 1:obj.dim
634 for l = 1:obj.dim
635 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
636 end
637 end
638
639 varargout{k} = alpha;
640
641 case 'alpha_tot'
642 % alpha = alpha(i,j) is the penalty strength for displacement BC.
643 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary);
644 E = obj.E;
645 [m, n] = size(alpha{1,1});
646 alpha_tot = sparse(m*obj.dim, n*obj.dim);
647 for i = 1:obj.dim
648 for l = 1:obj.dim
649 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
650 end
651 end
652 varargout{k} = alpha_tot;
653
501 otherwise 654 otherwise
502 error(['No such operator: operator = ' op{i}]); 655 error(['No such operator: operator = ' op{k}]);
503 end 656 end
504 end 657 end
505 658
659 end
660
661 function H = getBoundaryQuadrature(obj, boundary)
662 switch boundary
663 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
664 j = 1;
665 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
666 j = 2;
667 otherwise
668 error('No such boundary: boundary = %s',boundary);
669 end
670 H = obj.H_boundary{j};
671 I_dim = speye(obj.dim, obj.dim);
672 H = kron(H, I_dim);
506 end 673 end
507 674
508 function N = size(obj) 675 function N = size(obj)
509 N = obj.dim*prod(obj.m); 676 N = obj.dim*prod(obj.m);
510 end 677 end