Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 963:c75ddd568fcc feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Turn alpha into a boundary operator. Add properties H_w etc for getBoundaryQuadrature to work.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 19 Dec 2018 06:58:10 +0100
parents 2efeedf8c34f
children db3411264b96
comparison
equal deleted inserted replaced
962:262b52c3f268 963:c75ddd568fcc
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, H_1D % Inner products 33 H, Hi, H_1D % Inner products
34 e_l, e_r 34 e_l, e_r
35 e_w, e_e, e_s, e_n
36
37
35 d1_l, d1_r % Normal derivatives at the boundary 38 d1_l, d1_r % Normal derivatives at the boundary
36 E % E{i}^T picks out component i 39 E % E{i}^T picks out component i
37 40
38 H_boundary % Boundary inner products 41 H_boundary % Boundary inner products
42 H_w, H_e, H_s, H_n
39 43
40 % Kroneckered norms and coefficients 44 % Kroneckered norms and coefficients
41 RHOi_kron 45 RHOi_kron
42 Hi_kron 46 Hi_kron
43 47
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});
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};
186 obj.H_1D = {H{1}, H{2}}; 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};
187 199
188 % E{i}^T picks out component i. 200 % E{i}^T picks out component i.
189 E = cell(dim,1); 201 E = cell(dim,1);
190 I = speye(m_tot,m_tot); 202 I = speye(m_tot,m_tot);
191 for i = 1:dim 203 for i = 1:dim
328 340
329 % j is the coordinate direction of the boundary 341 % j is the coordinate direction of the boundary
330 j = obj.get_boundary_number(boundary); 342 j = obj.get_boundary_number(boundary);
331 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 343 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
332 344
345
333 E = obj.E; 346 E = obj.E;
334 Hi = obj.Hi; 347 Hi = obj.Hi;
335 LAMBDA = obj.LAMBDA; 348 LAMBDA = obj.LAMBDA;
336 MU = obj.MU; 349 MU = obj.MU;
337 RHOi = obj.RHOi; 350 RHOi = obj.RHOi;
352 alpha = obj.get_boundary_operator('alpha', boundary); 365 alpha = obj.get_boundary_operator('alpha', boundary);
353 366
354 % Loop over components that Dirichlet penalties end up on 367 % Loop over components that Dirichlet penalties end up on
355 for i = 1:dim 368 for i = 1:dim
356 C = transpose(T{k,i}); 369 C = transpose(T{k,i});
357 A = -alpha{i,k}; 370 A = -e*transpose(alpha{i,k});
358 B = A + e*C; 371 B = A + e*C;
359 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); 372 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
360 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 373 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
361 end 374 end
362 375
390 error('Unknown type of interpolation: %s ', type.interpolation); 403 error('Unknown type of interpolation: %s ', type.interpolation);
391 end 404 end
392 end 405 end
393 406
394 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 407 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
408 tuning = type.tuning;
409
395 % u denotes the solution in the own domain 410 % u denotes the solution in the own domain
396 % v denotes the solution in the neighbour domain 411 % v denotes the solution in the neighbour domain
397 % Operators without subscripts are from the own domain. 412 % Operators without subscripts are from the own domain.
398 413
399 % j is the coordinate direction of the boundary 414 % j is the coordinate direction of the boundary
455 % Loop over components that penalties end up on 470 % Loop over components that penalties end up on
456 for i = 1:dim 471 for i = 1:dim
457 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}'; 472 closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
458 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}'; 473 penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
459 474
460 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}; 475 closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau{i}';
461 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}; 476 penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau_v{i}';
462 477
463 % Loop over components that we have interface conditions on 478 % Loop over components that we have interface conditions on
464 for k = 1:dim 479 for k = 1:dim
465 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}'; 480 closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}';
466 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}'; 481 penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}';
523 case {'w','W','west','West','s','S','south','South'} 538 case {'w','W','west','West','s','S','south','South'}
524 varargout{k} = obj.d1_l{j}; 539 varargout{k} = obj.d1_l{j};
525 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 540 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
526 varargout{k} = obj.d1_r{j}; 541 varargout{k} = obj.d1_r{j};
527 end 542 end
528 case 'H'
529 varargout{k} = obj.H_boundary{j};
530 case 'T' 543 case 'T'
531 switch boundary 544 switch boundary
532 case {'w','W','west','West','s','S','south','South'} 545 case {'w','W','west','West','s','S','south','South'}
533 varargout{k} = obj.T_l{j}; 546 varargout{k} = obj.T_l{j};
534 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 547 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
539 case {'w','W','west','West','s','S','south','South'} 552 case {'w','W','west','West','s','S','south','South'}
540 varargout{k} = obj.tau_l{j}; 553 varargout{k} = obj.tau_l{j};
541 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 554 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
542 varargout{k} = obj.tau_r{j}; 555 varargout{k} = obj.tau_r{j};
543 end 556 end
557 case 'H'
558 varargout{k} = obj.H_boundary{j};
544 case 'alpha' 559 case 'alpha'
545 % alpha = alpha(i,j) is the penalty strength for displacement BC. 560 % alpha = alpha(i,j) is the penalty strength for displacement BC.
561 switch boundary
562 case {'w','W','west','West','s','S','south','South'}
563 e = obj.e_l{j};
564 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
565 e = obj.e_r{j};
566 end
567
546 tuning = 1.2; 568 tuning = 1.2;
547 LAMBDA = obj.LAMBDA; 569 LAMBDA = obj.LAMBDA;
548 MU = obj.MU; 570 MU = obj.MU;
549 571
550 dim = obj.dim; 572 dim = obj.dim;
563 alpha_func = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 585 alpha_func = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
564 + d(i,j)* a_mu_i*MU ... 586 + d(i,j)* a_mu_i*MU ...
565 + db(i,j)*a_mu_ij*MU ); 587 + db(i,j)*a_mu_ij*MU );
566 for i = 1:obj.dim 588 for i = 1:obj.dim
567 for l = 1:obj.dim 589 for l = 1:obj.dim
568 alpha{i,l} = d(i,l)*alpha_func(i,j); 590 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
569 end 591 end
570 end 592 end
571 593
572 varargout{k} = alpha; 594 varargout{k} = alpha;
573 otherwise 595 otherwise