Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 974:1c334842bf23 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Extract tuning from alpha.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 25 Dec 2018 16:39:06 +0100
parents f5a99cca4654
children b4fa176b4287
comparison
equal deleted inserted replaced
973:f5a99cca4654 974:1c334842bf23
382 alpha = obj.getBoundaryOperator('alpha', boundary); 382 alpha = obj.getBoundaryOperator('alpha', boundary);
383 383
384 % Loop over components that Dirichlet penalties end up on 384 % Loop over components that Dirichlet penalties end up on
385 for i = 1:dim 385 for i = 1:dim
386 C = transpose(T{k,i}); 386 C = transpose(T{k,i});
387 A = -e*transpose(alpha{i,k}); 387 A = -tuning*e*transpose(alpha{i,k});
388 B = A + e*C; 388 B = A + e*C;
389 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); 389 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
390 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 390 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
391 end 391 end
392 392
436 % Operators and quantities that correspond to the own domain only 436 % Operators and quantities that correspond to the own domain only
437 Hi = obj.Hi_kron; 437 Hi = obj.Hi_kron;
438 RHOi = obj.RHOi_kron; 438 RHOi = obj.RHOi_kron;
439 439
440 % Penalty strength operators 440 % Penalty strength operators
441 alpha_u = 1/4*obj.getBoundaryOperator('alpha_tot', boundary); 441 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary);
442 alpha_v = 1/4*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary); 442 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary);
443 443
444 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e'); 444 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
445 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v'); 445 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
446 446
447 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau'; 447 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
550 550
551 case 'alpha' 551 case 'alpha'
552 % alpha = alpha(i,j) is the penalty strength for displacement BC. 552 % alpha = alpha(i,j) is the penalty strength for displacement BC.
553 e = obj.getBoundaryOperator('e', boundary); 553 e = obj.getBoundaryOperator('e', boundary);
554 554
555 tuning = 1.2;
556 LAMBDA = obj.LAMBDA; 555 LAMBDA = obj.LAMBDA;
557 MU = obj.MU; 556 MU = obj.MU;
558 557
559 dim = obj.dim; 558 dim = obj.dim;
560 theta_R = obj.theta_R{j}; 559 theta_R = obj.theta_R{j};
567 566
568 d = @kroneckerDelta; % Kronecker delta 567 d = @kroneckerDelta; % Kronecker delta
569 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 568 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
570 alpha = cell(obj.dim, obj.dim); 569 alpha = cell(obj.dim, obj.dim);
571 570
572 alpha_func = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 571 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
573 + d(i,j)* a_mu_i*MU ... 572 + d(i,j)* a_mu_i*MU ...
574 + db(i,j)*a_mu_ij*MU ); 573 + db(i,j)*a_mu_ij*MU;
575 for i = 1:obj.dim 574 for i = 1:obj.dim
576 for l = 1:obj.dim 575 for l = 1:obj.dim
577 alpha{i,l} = d(i,l)*alpha_func(i,j)*e; 576 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
578 end 577 end
579 end 578 end