Mercurial > repos > public > sbplib

comparison +scheme/elasticShearVariable.m @ 673:9e1d2351f539 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Change penalties in Free BC to enable mixed Dirichlet/Free BC.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 23 Dec 2017 15:42:11 +0100
parents 17e62551cdc2
children eeaf9a00e304
comparison
equal deleted inserted replaced
670:52a9dedb9171 673:9e1d2351f539
6 grid 6 grid
7 dim 7 dim
8 8
9 order % Order accuracy for the approximation 9 order % Order accuracy for the approximation
10 10
11 A % Variable coefficient lambda of the operator (as diagonal matrix here) 11 A % Variable coefficient mu of the operator (as diagonal matrix here)
12 RHO % Density (as diagonal matrix here) 12 RHO % Density (as diagonal matrix here)
13 13
14 D % Total operator 14 D % Total operator
15 D1 % First derivatives 15 D1 % First derivatives
16 D2 % Second derivatives 16 D2 % Second derivatives
29 end 29 end
30 30
31 methods 31 methods
32 % Implements the shear part of the elastic wave equation, i.e. 32 % Implements the shear part of the elastic wave equation, i.e.
33 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i 33 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i
34 % where a = lambda. 34 % where a = mu.
35 35
36 function obj = elasticShearVariable(g ,order, a_fun, rho_fun, opSet) 36 function obj = elasticShearVariable(g ,order, a_fun, rho_fun, opSet)
37 default_arg('opSet',@sbp.D2Variable); 37 default_arg('opSet',@sbp.D2Variable);
38 default_arg('a_fun', @(x,y) 0*x+1); 38 default_arg('a_fun', @(x,y) 0*x+1);
39 default_arg('rho_fun', @(x,y) 0*x+1); 39 default_arg('rho_fun', @(x,y) 0*x+1);
184 end 184 end
185 185
186 186
187 % Closure functions return the operators applied to the own domain to close the boundary 187 % Closure functions return the operators applied to the own domain to close the boundary
188 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 188 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
189 % Here penalty{i,j} enforces data component j on solution component i
189 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 190 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
190 % type is a string specifying the type of boundary condition if there are several. 191 % type is a string specifying the type of boundary condition if there are several.
191 % data is a function returning the data that should be applied at the boundary. 192 % data is a function returning the data that should be applied at the boundary.
192 % neighbour_scheme is an instance of Scheme that should be interfaced to. 193 % neighbour_scheme is an instance of Scheme that should be interfaced to.
193 % neighbour_boundary is a string specifying which boundary to interface to. 194 % neighbour_boundary is a string specifying which boundary to interface to.
305 penalty = penalties; 306 penalty = penalties;
306 307
307 % Free boundary condition 308 % Free boundary condition
308 case {'F','f','Free','free'} 309 case {'F','f','Free','free'}
309 closures = cell(obj.dim,1); 310 closures = cell(obj.dim,1);
310 penalties = cell(obj.dim,1); 311 penalties = cell(obj.dim,obj.dim);
311 % Loop over components 312 % Loop over components
312 for i = 1:obj.dim 313 for i = 1:obj.dim
313 closures{i} = E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(... 314 closures{i} = E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(...
314 e{j}'*A*e{j}*d{j}'*E{i}' + ... 315 e{j}'*A*e{j}*d{j}'*E{i}' + ...
315 delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ... 316 delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ...
316 delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ... 317 delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ...
317 ); 318 );
318 penalties{i} = - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma; 319 penalties{i,i} = - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma;
319 end 320 end
320 [rows, cols] = size(closures{1}); 321 [rows, cols] = size(closures{1});
321 closure = sparse(rows, cols); 322 closure = sparse(rows, cols);
322 for i = 1:obj.dim 323 for i = 1:obj.dim
323 closure = closure + closures{i}; 324 closure = closure + closures{i};
325 for j = 1:obj.dim
326 if i~=j
327 [rows cols] = size(penalties{j,j});
328 penalties{i,j} = sparse(rows,cols);
329 end
330 end
324 end 331 end
325 penalty = penalties; 332 penalty = penalties;
326 333
327 334
328 % Unknown boundary condition 335 % Unknown boundary condition