Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 726:37d5d69b1a0d feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Refactor BC code in Elastic2dVariable
author Martin Almquist <malmquist@stanford.edu>
date Thu, 19 Apr 2018 17:06:27 -0700
parents 60eb7f46d8d9
children 6d5953fc090e
comparison
equal deleted inserted replaced
725:e9e15d64f9f9 726:37d5d69b1a0d
269 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 269 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
270 default_arg('type',{'free','free'}); 270 default_arg('type',{'free','free'});
271 default_arg('parameter', []); 271 default_arg('parameter', []);
272 272
273 % j is the coordinate direction of the boundary 273 % j is the coordinate direction of the boundary
274 % nj: outward unit normal component. 274 j = obj.get_boundary_number(boundary);
275 % nj = -1 for west, south, bottom boundaries 275 [e, T, tau] = obj.get_boundary_operator({'e','T','tau'}, boundary);
276 % nj = 1 for east, north, top boundaries
277 [j, nj] = obj.get_boundary_number(boundary);
278 switch nj
279 case 1
280 e = obj.e_r;
281 d = obj.d1_r;
282 tau = obj.tau_r{j};
283 T = obj.T_r{j};
284 case -1
285 e = obj.e_l;
286 d = obj.d1_l;
287 tau = obj.tau_l{j};
288 T = obj.T_l{j};
289 end
290 276
291 E = obj.E; 277 E = obj.E;
292 Hi = obj.Hi; 278 Hi = obj.Hi;
293 H_gamma = obj.H_boundary{j}; 279 H_gamma = obj.H_boundary{j};
294 LAMBDA = obj.LAMBDA; 280 LAMBDA = obj.LAMBDA;
334 % Loop over components that Dirichlet penalties end up on 320 % Loop over components that Dirichlet penalties end up on
335 for i = 1:dim 321 for i = 1:dim
336 C = T{k,i}; 322 C = T{k,i};
337 A = -d(i,k)*alpha(i,j); 323 A = -d(i,k)*alpha(i,j);
338 B = A + C; 324 B = A + C;
339 closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); 325 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
340 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; 326 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e*H_gamma;
341 end 327 end
342 328
343 % Free boundary condition 329 % Free boundary condition
344 case {'F','f','Free','free','traction','Traction','t','T'} 330 case {'F','f','Free','free','traction','Traction','t','T'}
345 closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); 331 closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} );
346 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; 332 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e*H_gamma;
347 333
348 % Unknown boundary condition 334 % Unknown boundary condition
349 otherwise 335 otherwise
350 error('No such boundary condition: type = %s',type); 336 error('No such boundary condition: type = %s',type);
351 end 337 end
378 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 364 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
379 nj = 1; 365 nj = 1;
380 end 366 end
381 end 367 end
382 368
383 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 369 % Returns the boundary operator op for the boundary specified by the string boundary.
384 function [return_op] = get_boundary_operator(obj, op, boundary) 370 % op: may be a cell array of strings
371 function [varargout] = get_boundary_operator(obj, op, boundary)
385 372
386 switch boundary 373 switch boundary
387 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 374 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
388 j = 1; 375 j = 1;
389 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 376 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
390 j = 2; 377 j = 2;
391 otherwise 378 otherwise
392 error('No such boundary: boundary = %s',boundary); 379 error('No such boundary: boundary = %s',boundary);
393 end 380 end
394 381
395 switch op 382 if ~iscell(op)
396 case 'e' 383 op = {op};
397 switch boundary 384 end
398 case {'w','W','west','West','s','S','south','South'} 385
399 return_op = obj.e_l{j}; 386 for i = 1:length(op)
400 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 387 switch op{i}
401 return_op = obj.e_r{j}; 388 case 'e'
402 end 389 switch boundary
403 case 'd' 390 case {'w','W','west','West','s','S','south','South'}
404 switch boundary 391 varargout{i} = obj.e_l{j};
405 case {'w','W','west','West','s','S','south','South'} 392 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
406 return_op = obj.d1_l{j}; 393 varargout{i} = obj.e_r{j};
407 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 394 end
408 return_op = obj.d1_r{j}; 395 case 'd'
409 end 396 switch boundary
410 otherwise 397 case {'w','W','west','West','s','S','south','South'}
411 error(['No such operator: operatr = ' op]); 398 varargout{i} = obj.d1_l{j};
399 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
400 varargout{i} = obj.d1_r{j};
401 end
402 case 'H'
403 varargout{i} = obj.H_boundary{j};
404 case 'T'
405 switch boundary
406 case {'w','W','west','West','s','S','south','South'}
407 varargout{i} = obj.T_l{j};
408 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
409 varargout{i} = obj.T_r{j};
410 end
411 case 'tau'
412 switch boundary
413 case {'w','W','west','West','s','S','south','South'}
414 varargout{i} = obj.tau_l{j};
415 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
416 varargout{i} = obj.tau_r{j};
417 end
418 otherwise
419 error(['No such operator: operator = ' op{i}]);
420 end
412 end 421 end
413 422
414 end 423 end
415 424
416 function N = size(obj) 425 function N = size(obj)