Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.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 aa4ef495f1fd
children b374a8aa9246 5751262b323b
comparison
equal deleted inserted replaced
794:0cac17097c37 795:1f6b2fb69225
267 267
268 268
269 % Closure functions return the operators applied to the own domain to close the boundary 269 % Closure functions return the operators applied to the own domain to close the boundary
270 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 270 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
271 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 271 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
272 % type is a cell array of strings specifying the type of boundary condition for each component. 272 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
273 % on the first component.
273 % data is a function returning the data that should be applied at the boundary. 274 % data is a function returning the data that should be applied at the boundary.
274 % neighbour_scheme is an instance of Scheme that should be interfaced to. 275 % neighbour_scheme is an instance of Scheme that should be interfaced to.
275 % neighbour_boundary is a string specifying which boundary to interface to. 276 % neighbour_boundary is a string specifying which boundary to interface to.
276 function [closure, penalty] = boundary_condition(obj, boundary, type, tuning) 277 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
277 default_arg('type',{'free','free'});
278 default_arg('tuning', 1.2); 278 default_arg('tuning', 1.2);
279 279
280 if ~iscell(type) 280 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
281 type = {type, type}; 281 comp = bc{1};
282 end 282 type = bc{2};
283 283
284 % j is the coordinate direction of the boundary 284 % j is the coordinate direction of the boundary
285 j = obj.get_boundary_number(boundary); 285 j = obj.get_boundary_number(boundary);
286 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); 286 [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
287 287
294 dim = obj.dim; 294 dim = obj.dim;
295 m_tot = obj.grid.N(); 295 m_tot = obj.grid.N();
296 296
297 % Preallocate 297 % Preallocate
298 closure = sparse(dim*m_tot, dim*m_tot); 298 closure = sparse(dim*m_tot, dim*m_tot);
299 penalty = cell(dim,1); 299 penalty = sparse(dim*m_tot, m_tot/obj.m(j));
300 for k = 1:dim 300
301 penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j)); 301 k = comp;
302 end 302 switch type
303 303
304 % Loop over components that we (potentially) have different BC on 304 % Dirichlet boundary condition
305 for k = 1:dim 305 case {'D','d','dirichlet','Dirichlet'}
306 switch type{k} 306
307 307 phi = obj.phi{j};
308 % Dirichlet boundary condition 308 h = obj.h(j);
309 case {'D','d','dirichlet','Dirichlet'} 309 h11 = obj.H11{j}*h;
310 310 gamma = obj.gamma{j};
311 phi = obj.phi{j}; 311
312 h = obj.h(j); 312 a_lambda = dim/h11 + 1/(h11*phi);
313 h11 = obj.H11{j}*h; 313 a_mu_i = 2/(gamma*h);
314 gamma = obj.gamma{j}; 314 a_mu_ij = 2/h11 + 1/(h11*phi);
315 315
316 a_lambda = dim/h11 + 1/(h11*phi); 316 d = @kroneckerDelta; % Kronecker delta
317 a_mu_i = 2/(gamma*h); 317 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
318 a_mu_ij = 2/h11 + 1/(h11*phi); 318 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
319 319 + d(i,j)* a_mu_i*MU ...
320 d = @kroneckerDelta; % Kronecker delta 320 + db(i,j)*a_mu_ij*MU );
321 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 321
322 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 322 % Loop over components that Dirichlet penalties end up on
323 + d(i,j)* a_mu_i*MU ... 323 for i = 1:dim
324 + db(i,j)*a_mu_ij*MU ); 324 C = T{k,i};
325 325 A = -d(i,k)*alpha(i,j);
326 % Loop over components that Dirichlet penalties end up on 326 B = A + C;
327 for i = 1:dim 327 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
328 C = T{k,i}; 328 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
329 A = -d(i,k)*alpha(i,j); 329 end
330 B = A + C; 330
331 closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); 331 % Free boundary condition
332 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e*H_gamma; 332 case {'F','f','Free','free','traction','Traction','t','T'}
333 end 333 closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} );
334 334 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
335 % Free boundary condition 335
336 case {'F','f','Free','free','traction','Traction','t','T'} 336 % Unknown boundary condition
337 closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); 337 otherwise
338 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e*H_gamma; 338 error('No such boundary condition: type = %s',type);
339
340 % Unknown boundary condition
341 otherwise
342 error('No such boundary condition: type = %s',type);
343 end
344 end 339 end
345 end 340 end
346 341
347 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 342 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
348 % u denotes the solution in the own domain 343 % u denotes the solution in the own domain