Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1203:25cadc69a589 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update comments in Anisotropic
author Martin Almquist <malmquist@stanford.edu>
date Thu, 05 Sep 2019 17:10:58 -0700
parents 31d7288d0653
children 687515778437
comparison
equal deleted inserted replaced
1202:31d7288d0653 1203:25cadc69a589
365 % {'normal', 'd'} or {'tangential', 't'} for conditions on 365 % {'normal', 'd'} or {'tangential', 't'} for conditions on
366 % tangential/normal component. 366 % tangential/normal component.
367 % data is a function returning the data that should be applied at the boundary. 367 % data is a function returning the data that should be applied at the boundary.
368 % neighbour_scheme is an instance of Scheme that should be interfaced to. 368 % neighbour_scheme is an instance of Scheme that should be interfaced to.
369 % neighbour_boundary is a string specifying which boundary to interface to. 369 % neighbour_boundary is a string specifying which boundary to interface to.
370
371 % For displacement bc:
372 % bc = {comp, 'd', dComps},
373 % where
374 % dComps = vector of components with displacement BC. Default: 1:dim.
375 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
370 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) 376 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
371 default_arg('tuning', 1.0); 377 default_arg('tuning', 1.0);
372 378
373 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); 379 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
374 comp = bc{1}; 380 comp = bc{1};
414 % Y: symmetrizing part of penalty 420 % Y: symmetrizing part of penalty
415 % Z: symmetric part of penalty 421 % Z: symmetric part of penalty
416 % X = Y + Z. 422 % X = Y + Z.
417 423
418 % Nonsymmetric part goes on all components to 424 % Nonsymmetric part goes on all components to
419 % Yield traction in discrete energy rate 425 % yield traction in discrete energy rate
420 for i = 1:dim 426 for i = 1:dim
421 Y = T{j,i}'; 427 Y = T{j,i}';
422 X = e*Y; 428 X = e*Y;
423 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); 429 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
424 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; 430 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
425 end 431 end
426 432
427 % Symmetric part only required on components with displacement BC 433 % Symmetric part only required on components with displacement BC.
434 % (Otherwise it's not symmetric.)
428 for i = dComps 435 for i = dComps
429 Z = sparse(m_tot, m_tot); 436 Z = sparse(m_tot, m_tot);
430 for l = 1:dim 437 for l = 1:dim
431 for k = 1:dim 438 for k = 1:dim
432 Z = Z + nu(l)*C{i,l,k,j}*nu(k); 439 Z = Z + nu(l)*C{i,l,k,j}*nu(k);