Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 905:459eeb99130f feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Include type as (optional) input parameter in the interface method of all schemes.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 22 Nov 2018 22:03:44 -0800
parents 60eb7f46d8d9
children b9c98661ff5d
comparison
equal deleted inserted replaced
904:14b093a344eb 905:459eeb99130f
1 classdef Elastic2dVariable < scheme.Scheme 1 classdef Elastic2dVariable < scheme.Scheme
2 2
3 % Discretizes the elastic wave equation: 3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i 4 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
5 % opSet should be cell array of opSets, one per dimension. This 5 % opSet should be cell array of opSets, one per dimension. This
6 % is useful if we have periodic BC in one direction. 6 % is useful if we have periodic BC in one direction.
7 7
8 properties 8 properties
9 m % Number of points in each direction, possibly a vector 9 m % Number of points in each direction, possibly a vector
35 gamma % Borrowing constant for d1 from M 35 gamma % Borrowing constant for d1 from M
36 H11 % First element of H 36 H11 % First element of H
37 e_l, e_r 37 e_l, e_r
38 d1_l, d1_r % Normal derivatives at the boundary 38 d1_l, d1_r % Normal derivatives at the boundary
39 E % E{i}^T picks out component i 39 E % E{i}^T picks out component i
40 40
41 H_boundary % Boundary inner products 41 H_boundary % Boundary inner products
42 42
43 % Kroneckered norms and coefficients 43 % Kroneckered norms and coefficients
44 RHOi_kron 44 RHOi_kron
45 Hi_kron 45 Hi_kron
221 % Loop over components 221 % Loop over components
222 for i = 1:dim 222 for i = 1:dim
223 tau_l{j}{i} = sparse(m_tot,dim*m_tot); 223 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
224 tau_r{j}{i} = sparse(m_tot,dim*m_tot); 224 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
225 for k = 1:dim 225 for k = 1:dim
226 T_l{j}{i,k} = ... 226 T_l{j}{i,k} = ...
227 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... 227 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
228 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... 228 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
229 -d(i,k)*MU*e_l{j}*d1_l{j}'; 229 -d(i,k)*MU*e_l{j}*d1_l{j}';
230 230
231 T_r{j}{i,k} = ... 231 T_r{j}{i,k} = ...
232 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... 232 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
233 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... 233 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
234 +d(i,k)*MU*e_r{j}*d1_r{j}'; 234 +d(i,k)*MU*e_r{j}*d1_r{j}';
235 235
236 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; 236 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
237 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; 237 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
238 end 238 end
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 % nj: outward unit normal component.
275 % nj = -1 for west, south, bottom boundaries 275 % nj = -1 for west, south, bottom boundaries
276 % nj = 1 for east, north, top boundaries 276 % nj = 1 for east, north, top boundaries
277 [j, nj] = obj.get_boundary_number(boundary); 277 [j, nj] = obj.get_boundary_number(boundary);
278 switch nj 278 switch nj
279 case 1 279 case 1
327 327
328 d = @kroneckerDelta; % Kronecker delta 328 d = @kroneckerDelta; % Kronecker delta
329 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 329 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
330 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... 330 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
331 + d(i,j)* a_mu_i*MU ... 331 + d(i,j)* a_mu_i*MU ...
332 + db(i,j)*a_mu_ij*MU ); 332 + db(i,j)*a_mu_ij*MU );
333 333
334 % Loop over components that Dirichlet penalties end up on 334 % Loop over components that Dirichlet penalties end up on
335 for i = 1:dim 335 for i = 1:dim
336 C = T{k,i}; 336 C = T{k,i};
337 A = -d(i,k)*alpha(i,j); 337 A = -d(i,k)*alpha(i,j);
338 B = A + C; 338 B = A + C;
339 closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); 339 closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' );
340 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; 340 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma;
341 end 341 end
342 342
343 % Free boundary condition 343 % Free boundary condition
344 case {'F','f','Free','free','traction','Traction','t','T'} 344 case {'F','f','Free','free','traction','Traction','t','T'}
345 closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} ); 345 closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} );
346 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; 346 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma;
347 347
348 % Unknown boundary condition 348 % Unknown boundary condition
349 otherwise 349 otherwise
350 error('No such boundary condition: type = %s',type); 350 error('No such boundary condition: type = %s',type);
351 end 351 end
352 end 352 end
353 end 353 end
354 354
355 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 355 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
356 % u denotes the solution in the own domain 356 % u denotes the solution in the own domain
357 % v denotes the solution in the neighbour domain 357 % v denotes the solution in the neighbour domain
358 tuning = 1.2; 358 tuning = 1.2;
359 % tuning = 20.2; 359 % tuning = 20.2;
360 error('Interface not implemented'); 360 error('Interface not implemented');