Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinear.m @ 795:1f6b2fb69225 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Revert bcSetup and update bc functions in elastic schemes to be compatible.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 25 Jul 2018 18:33:03 -0700
parents 08f3ffe63f48
children cb4bfd0d81ea b8bd54332763
comparison
equal deleted inserted replaced
794:0cac17097c37 795:1f6b2fb69225
371 371
372 372
373 % Closure functions return the operators applied to the own domain to close the boundary 373 % Closure functions return the operators applied to the own domain to close the boundary
374 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 374 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
375 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 375 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
376 % type is a cell array of strings specifying the type of boundary condition for each component. 376 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
377 % on the first component.
377 % data is a function returning the data that should be applied at the boundary. 378 % data is a function returning the data that should be applied at the boundary.
378 % neighbour_scheme is an instance of Scheme that should be interfaced to. 379 % neighbour_scheme is an instance of Scheme that should be interfaced to.
379 % neighbour_boundary is a string specifying which boundary to interface to. 380 % neighbour_boundary is a string specifying which boundary to interface to.
380 function [closure, penalty] = boundary_condition(obj, boundary, type, tuning) 381 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
381 default_arg('type',{'free','free'});
382 default_arg('tuning', 1.2); 382 default_arg('tuning', 1.2);
383 383
384 if ~iscell(type) 384 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
385 type = {type, type}; 385 comp = bc{1};
386 end 386 type = bc{2};
387 387
388 % j is the coordinate direction of the boundary 388 % j is the coordinate direction of the boundary
389 j = obj.get_boundary_number(boundary); 389 j = obj.get_boundary_number(boundary);
390 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 390 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
391 391
399 dim = obj.dim; 399 dim = obj.dim;
400 m_tot = obj.grid.N(); 400 m_tot = obj.grid.N();
401 401
402 % Preallocate 402 % Preallocate
403 closure = sparse(dim*m_tot, dim*m_tot); 403 closure = sparse(dim*m_tot, dim*m_tot);
404 penalty = cell(dim,1); 404 penalty = sparse(dim*m_tot, m_tot/obj.m(j));
405 for k = 1:dim
406 penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
407 end
408 405
409 % Loop over components that we (potentially) have different BC on 406 % Loop over components that we (potentially) have different BC on
410 for k = 1:dim 407 k = comp;
411 switch type{k} 408 switch type
412 409
413 % Dirichlet boundary condition 410 % Dirichlet boundary condition
414 case {'D','d','dirichlet','Dirichlet'} 411 case {'D','d','dirichlet','Dirichlet'}
415 412
416 phi = obj.phi{j}; 413 phi = obj.phi{j};
417 h = obj.h(j); 414 h = obj.h(j);
418 h11 = obj.H11{j}*h; 415 h11 = obj.H11{j}*h;
419 gamma = obj.gamma{j}; 416 gamma = obj.gamma{j};
420 417
421 a_lambda = dim/h11 + 1/(h11*phi); 418 a_lambda = dim/h11 + 1/(h11*phi);
422 a_mu_i = 2/(gamma*h); 419 a_mu_i = 2/(gamma*h);
423 a_mu_ij = 2/h11 + 1/(h11*phi); 420 a_mu_ij = 2/h11 + 1/(h11*phi);
424 421
425 d = @kroneckerDelta; % Kronecker delta 422 d = @kroneckerDelta; % Kronecker delta
426 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 423 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
427 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 424 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
428 + d(i,j)* a_mu_i*MU ... 425 + d(i,j)* a_mu_i*MU ...
429 + db(i,j)*a_mu_ij*MU ); 426 + db(i,j)*a_mu_ij*MU );
430 427
431 % Loop over components that Dirichlet penalties end up on 428 % Loop over components that Dirichlet penalties end up on
432 for i = 1:dim 429 for i = 1:dim
433 C = T{k,i}; 430 C = T{k,i};
434 A = -d(i,k)*alpha(i,j); 431 A = -d(i,k)*alpha(i,j);
435 B = A + C; 432 B = A + C;
436 closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' ); 433 closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' );
437 penalty{k} = penalty{k} - E{i}*RHOi*Hi*Ji*B'*e*H_gamma; 434 penalty = penalty - E{i}*RHOi*Hi*Ji*B'*e*H_gamma;
438 end 435 end
439 436
440 % Free boundary condition 437 % Free boundary condition
441 case {'F','f','Free','free','traction','Traction','t','T'} 438 case {'F','f','Free','free','traction','Traction','t','T'}
442 closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} ); 439 closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} );
443 penalty{k} = penalty{k} + E{k}*RHOi*Ji*Hi*e*H_gamma; 440 penalty = penalty + E{k}*RHOi*Ji*Hi*e*H_gamma;
444 441
445 % Unknown boundary condition 442 % Unknown boundary condition
446 otherwise 443 otherwise
447 error('No such boundary condition: type = %s',type); 444 error('No such boundary condition: type = %s',type);
448 end
449 end 445 end
450 end 446 end
451 447
452 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 448 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
453 % u denotes the solution in the own domain 449 % u denotes the solution in the own domain