Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariableAnisotropic.m @ 1202:31d7288d0653 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make Anisotropic work for mixed displacement/traction BC at the same boundary.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 05 Sep 2019 16:41:49 -0700
parents 8f4e79aa32ba
children 25cadc69a589
comparison
equal deleted inserted replaced
1201:8f4e79aa32ba 1202:31d7288d0653
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 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) 370 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
371 default_arg('tuning', 1.0); 371 default_arg('tuning', 1.0);
372 372
373 assert( iscell(bc), 'The BC type must be a 2x1 cell array' ); 373 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
374 comp = bc{1}; 374 comp = bc{1};
375 type = bc{2}; 375 type = bc{2};
376 if ischar(comp) 376 if ischar(comp)
377 comp = obj.getComponent(comp, boundary); 377 comp = obj.getComponent(comp, boundary);
378 end 378 end
402 402
403 % Dirichlet boundary condition 403 % Dirichlet boundary condition
404 % OBS! Cannot yet set one component at a time unless one assumes Displacement for all components 404 % OBS! Cannot yet set one component at a time unless one assumes Displacement for all components
405 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} 405 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
406 406
407 % Loop over components that Dirichlet penalties end up on 407 if numel(bc) >= 3
408 dComps = bc{3};
409 else
410 dComps = 1:dim;
411 end
412
413 % Loops over components that Dirichlet penalties end up on
408 % Y: symmetrizing part of penalty 414 % Y: symmetrizing part of penalty
409 % Z: symmetric part of penalty 415 % Z: symmetric part of penalty
410 % X = Y + Z. 416 % X = Y + Z.
417
418 % Nonsymmetric part goes on all components to
419 % Yield traction in discrete energy rate
411 for i = 1:dim 420 for i = 1:dim
412 Y = T{j,i}'; 421 Y = T{j,i}';
413 422 X = e*Y;
423 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
424 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
425 end
426
427 % Symmetric part only required on components with displacement BC
428 for i = dComps
414 Z = sparse(m_tot, m_tot); 429 Z = sparse(m_tot, m_tot);
415 for l = 1:dim 430 for l = 1:dim
416 for k = 1:dim 431 for k = 1:dim
417 Z = Z + nu(l)*C{i,l,k,j}*nu(k); 432 Z = Z + nu(l)*C{i,l,k,j}*nu(k);
418 end 433 end
419 end 434 end
420 Z = -tuning*dim/h11*Z; 435 Z = -tuning*dim/h11*Z;
421 X = Z + e*Y; 436 X = Z;
422 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); 437 closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' );
423 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; 438 penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma;
424 end 439 end
425 440
426 % Free boundary condition 441 % Free boundary condition