Mercurial > repos > public > sbplib

comparison +scheme/elasticVariable.m @ 679:247b58a4dbe8 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add support for Dirichlet and Traction BC on different components at the same boundary. Remove some unused variables and improve comments.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 05 Feb 2018 14:45:26 -0800
parents 06676c40e77f
children b035902869a8
comparison
equal deleted inserted replaced
678:06676c40e77f 679:247b58a4dbe8
1 classdef elasticVariable < scheme.Scheme 1 classdef elasticVariable < 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
6 5
7 properties 6 properties
8 m % Number of points in each direction, possibly a vector 7 m % Number of points in each direction, possibly a vector
9 h % Grid spacing 8 h % Grid spacing
10 9
14 order % Order of accuracy for the approximation 13 order % Order of accuracy for the approximation
15 14
16 % Diagonal matrices for varible coefficients 15 % Diagonal matrices for varible coefficients
17 LAMBDA % Variable coefficient, related to dilation 16 LAMBDA % Variable coefficient, related to dilation
18 MU % Shear modulus, variable coefficient 17 MU % Shear modulus, variable coefficient
19 RHO, RHOi % Density 18 RHO, RHOi % Density, variable
20 19
21 D % Total operator 20 D % Total operator
22 D1 % First derivatives 21 D1 % First derivatives
23 22
24 % Second derivatives 23 % Second derivatives
29 T_l, T_r 28 T_l, T_r
30 tau_l, tau_r 29 tau_l, tau_r
31 30
32 H, Hi % Inner products 31 H, Hi % Inner products
33 phi % Borrowing constant for (d1 - e^T*D1) from R 32 phi % Borrowing constant for (d1 - e^T*D1) from R
33 gamma % Borrowing constant for d1 from M
34 H11 % First element of H 34 H11 % First element of H
35 e_l, e_r 35 e_l, e_r
36 d1_l, d1_r % Normal derivatives at the boundary 36 d1_l, d1_r % Normal derivatives at the boundary
37 E % E{i}^T picks out component i 37 E % E{i}^T picks out component i
38 38
39 H_boundary % Boundary inner products 39 H_boundary % Boundary inner products
40 40
41 LAMBDA_boundary_l % Variable coefficients at boundaries
42 LAMBDA_boundary_r
43 MU_boundary_l
44 MU_boundary_r
45
46 % Kroneckered norms and coefficients 41 % Kroneckered norms and coefficients
47 RHOi_kron 42 RHOi_kron
48 Hi_kron 43 Hi_kron
49 end 44 end
50 45
51 methods 46 methods
52 % Implements the shear part of the elastic wave equation, i.e.
53 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i
54 % where a = mu.
55 47
56 function obj = elasticVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) 48 function obj = elasticVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
57 default_arg('opSet',@sbp.D2Variable); 49 default_arg('opSet',@sbp.D2Variable);
58 default_arg('lambda_fun', @(x,y) 0*x+1); 50 default_arg('lambda_fun', @(x,y) 0*x+1);
59 default_arg('mu_fun', @(x,y) 0*x+1); 51 default_arg('mu_fun', @(x,y) 0*x+1);
79 71
80 % Borrowing constants 72 % Borrowing constants
81 beta = ops{1}.borrowing.R.delta_D; 73 beta = ops{1}.borrowing.R.delta_D;
82 obj.H11 = ops{1}.borrowing.H11; 74 obj.H11 = ops{1}.borrowing.H11;
83 obj.phi = beta/obj.H11; 75 obj.phi = beta/obj.H11;
76 obj.gamma = ops{1}.borrowing.M.d1;
84 77
85 I = cell(dim,1); 78 I = cell(dim,1);
86 D1 = cell(dim,1); 79 D1 = cell(dim,1);
87 D2 = cell(dim,1); 80 D2 = cell(dim,1);
88 H = cell(dim,1); 81 H = cell(dim,1);
165 obj.H = kron(H{1},H{2}); 158 obj.H = kron(H{1},H{2});
166 obj.Hi = inv(obj.H); 159 obj.Hi = inv(obj.H);
167 obj.H_boundary = cell(dim,1); 160 obj.H_boundary = cell(dim,1);
168 obj.H_boundary{1} = H{2}; 161 obj.H_boundary{1} = H{2};
169 obj.H_boundary{2} = H{1}; 162 obj.H_boundary{2} = H{1};
170
171 % Boundary coefficient matrices
172 obj.LAMBDA_boundary_l = cell(dim,1);
173 obj.LAMBDA_boundary_r = cell(dim,1);
174 obj.MU_boundary_l = cell(dim,1);
175 obj.MU_boundary_r = cell(dim,1);
176 for i = 1:dim
177 obj.LAMBDA_boundary_l{i} = obj.e_l{i}'*LAMBDA*obj.e_l{i};
178 obj.LAMBDA_boundary_r{i} = obj.e_r{i}'*LAMBDA*obj.e_r{i};
179 obj.MU_boundary_l{i} = obj.e_l{i}'*MU*obj.e_l{i};
180 obj.MU_boundary_r{i} = obj.e_r{i}'*MU*obj.e_r{i};
181 end
182 163
183 % E{i}^T picks out component i. 164 % E{i}^T picks out component i.
184 E = cell(dim,1); 165 E = cell(dim,1);
185 I = speye(m_tot,m_tot); 166 I = speye(m_tot,m_tot);
186 for i = 1:dim 167 for i = 1:dim
275 256
276 257
277 % Closure functions return the operators applied to the own domain to close the boundary 258 % Closure functions return the operators applied to the own domain to close the boundary
278 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 259 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
279 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 260 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
280 % type is a string specifying the type of boundary condition if there are several. 261 % type is a cell array of strings specifying the type of boundary condition for each component.
281 % data is a function returning the data that should be applied at the boundary. 262 % data is a function returning the data that should be applied at the boundary.
282 % neighbour_scheme is an instance of Scheme that should be interfaced to. 263 % neighbour_scheme is an instance of Scheme that should be interfaced to.
283 % neighbour_boundary is a string specifying which boundary to interface to. 264 % neighbour_boundary is a string specifying which boundary to interface to.
284 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 265 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
285 default_arg('type','free'); 266 default_arg('type',{'free','free'});
286 default_arg('parameter', []); 267 default_arg('parameter', []);
287 268
288 % j is the coordinate direction of the boundary 269 % j is the coordinate direction of the boundary
289 % nj: outward unit normal component. 270 % nj: outward unit normal component.
290 % nj = -1 for west, south, bottom boundaries 271 % nj = -1 for west, south, bottom boundaries
308 H_gamma = obj.H_boundary{j}; 289 H_gamma = obj.H_boundary{j};
309 LAMBDA = obj.LAMBDA; 290 LAMBDA = obj.LAMBDA;
310 MU = obj.MU; 291 MU = obj.MU;
311 RHOi = obj.RHOi; 292 RHOi = obj.RHOi;
312 293
313 phi = obj.phi;
314 H11 = obj.H11;
315 h = obj.h;
316 dim = obj.dim; 294 dim = obj.dim;
317 m_tot = obj.grid.N(); 295 m_tot = obj.grid.N();
318 296
319 RHOi_kron = obj.RHOi_kron; 297 RHOi_kron = obj.RHOi_kron;
320 Hi_kron = obj.Hi_kron; 298 Hi_kron = obj.Hi_kron;
321 299
300 % Preallocate
322 closure = sparse(dim*m_tot, dim*m_tot); 301 closure = sparse(dim*m_tot, dim*m_tot);
323 penalty = cell(dim,1); 302 penalty = cell(dim,1);
324 switch type 303 for k = 1:dim
304 penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
305 end
306
307 % Loop over components that we (potentially) have different BC on
308 for k = 1:dim
309 switch type{k}
310
325 % Dirichlet boundary condition 311 % Dirichlet boundary condition
326 case {'D','d','dirichlet'} 312 case {'D','d','dirichlet','Dirichlet'}
327 error('not implemented') 313
328 tuning = 1.2; 314 tuning = 1.2;
329 phi = obj.phi; 315 phi = obj.phi;
330 316 h = obj.h(j);
331 sigma = tuning * obj.dim/(H11*h(j)) +... 317 h11 = obj.H11*h;
332 tuning * 1/(H11*h(j)*phi); 318 gamma = obj.gamma;
333 319
334 closure = - sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma*e{j}'*E{j}' ... 320 a_lambda = dim/h11 + 1/(h11*phi);
335 + nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma*e{j}'*E{j}'; 321 a_mu_i = 2/(gamma*h);
336 322 a_mu_ij = 2/h11 + 1/(h11*phi);
337 penalty = + sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma ... 323
338 - nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma; 324 d = @kroneckerDelta; % Kronecker delta
325 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
326 alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
327 + d(i,j)* a_mu_i*MU ...
328 + db(i,j)*a_mu_ij*MU );
329
330 % Loop over components that Dirichlet penalties end up on
331 for i = 1:dim
332 C = T{k,i};
333 A = -d(i,k)*alpha(i,j);
334 B = A + C;
335 closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' );
336 penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma;
337 end
339 338
340 % Free boundary condition 339 % Free boundary condition
341 case {'F','f','Free','free'} 340 case {'F','f','Free','free','traction','Traction','t','T'}
342 341 closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} );
343 % Loop over components of traction 342 penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma;
344 for i = 1:dim
345 closure = closure - E{i}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{i} );
346 penalty{i} = E{i}*RHOi*Hi*e{j}*H_gamma;
347 end
348
349 343
350 % Unknown boundary condition 344 % Unknown boundary condition
351 otherwise 345 otherwise
352 error('No such boundary condition: type = %s',type); 346 error('No such boundary condition: type = %s',type);
347 end
353 end 348 end
354 end 349 end
355 350
356 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) 351 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
357 % u denotes the solution in the own domain 352 % u denotes the solution in the own domain