Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 961:2efeedf8c34f feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make ElasticVariable2d.boundary_condition use get_boundary_operator for penalty strength.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 17 Dec 2018 21:07:59 -0800
parents ac566f3dc9b3
children c75ddd568fcc
comparison
equal deleted inserted replaced
960:ac566f3dc9b3 961:2efeedf8c34f
347 switch type 347 switch type
348 348
349 % Dirichlet boundary condition 349 % Dirichlet boundary condition
350 case {'D','d','dirichlet','Dirichlet'} 350 case {'D','d','dirichlet','Dirichlet'}
351 351
352 theta_R = obj.theta_R{j}; 352 alpha = obj.get_boundary_operator('alpha', boundary);
353 theta_H = obj.theta_H{j};
354 theta_M = obj.theta_M{j};
355
356 a_lambda = dim/theta_H + 1/theta_R;
357 a_mu_i = 2/theta_M;
358 a_mu_ij = 2/theta_H + 1/theta_R;
359
360 d = @kroneckerDelta; % Kronecker delta
361 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
362 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
363 + d(i,j)* a_mu_i*MU ...
364 + db(i,j)*a_mu_ij*MU );
365 353
366 % Loop over components that Dirichlet penalties end up on 354 % Loop over components that Dirichlet penalties end up on
367 for i = 1:dim 355 for i = 1:dim
368 C = transpose(T{k,i}); 356 C = transpose(T{k,i});
369 A = -d(i,k)*alpha(i,j); 357 A = -alpha{i,k};
370 B = A + e*C; 358 B = A + e*C;
371 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); 359 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
372 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 360 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
373 end 361 end
374 362