Mercurial > repos > public > sbplib

comparison +scheme/elasticDilationVariable.m @ 676:9926efb39330 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Clean up dilation BC code.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 18 Jan 2018 14:36:59 -0800
parents 90bf651abc7c
children
comparison
equal deleted inserted replaced
675:90bf651abc7c 676:9926efb39330
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 lambda of the operator (as diagonal matrix here)
12 RHO % Density (as diagonal matrix here) 12 RHO, RHOi % 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
17 Div % Divergence operator used for BC 17 Div % Divergence operator used for BC
24 24
25 H_boundary % Boundary inner products 25 H_boundary % Boundary inner products
26 26
27 A_boundary_l % Variable coefficient at boundaries 27 A_boundary_l % Variable coefficient at boundaries
28 A_boundary_r % 28 A_boundary_r %
29
30 % Kroneckered norms and coefficients
31 RHOi_kron
32 Hi_kron
29 end 33 end
30 34
31 methods 35 methods
32 % Implements the shear part of the elastic wave equation, i.e. 36 % 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 37 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i
85 %====== Assemble full operators ======== 89 %====== Assemble full operators ========
86 A = spdiag(a); 90 A = spdiag(a);
87 obj.A = A; 91 obj.A = A;
88 RHO = spdiag(rho); 92 RHO = spdiag(rho);
89 obj.RHO = RHO; 93 obj.RHO = RHO;
90 94 obj.RHOi = inv(RHO);
91 95
92 obj.D1 = cell(dim,1); 96 obj.D1 = cell(dim,1);
93 obj.D2 = cell(dim,1); 97 obj.D2 = cell(dim,1);
94 obj.e_l = cell(dim,1); 98 obj.e_l = cell(dim,1);
95 obj.e_r = cell(dim,1); 99 obj.e_r = cell(dim,1);
181 + db(i,j)*obj.D1{j}*E{j}'; 185 + db(i,j)*obj.D1{j}*E{j}';
182 end 186 end
183 end 187 end
184 obj.Div = Div; 188 obj.Div = Div;
185 189
190 % Kroneckered norms and coefficients
191 I_dim = speye(dim);
192 obj.RHOi_kron = kron(obj.RHOi, I_dim);
193 obj.Hi_kron = kron(obj.Hi, I_dim);
194
186 % Misc. 195 % Misc.
187 obj.m = m; 196 obj.m = m;
188 obj.h = h; 197 obj.h = h;
189 obj.order = order; 198 obj.order = order;
190 obj.grid = g; 199 obj.grid = g;
193 end 202 end
194 203
195 204
196 % Closure functions return the operators applied to the own domain to close the boundary 205 % Closure functions return the operators applied to the own domain to close the boundary
197 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 206 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
198 % Here penalty{i,j} enforces data component j on solution component i
199 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 207 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
200 % type is a string specifying the type of boundary condition if there are several. 208 % type is a string specifying the type of boundary condition if there are several.
201 % data is a function returning the data that should be applied at the boundary. 209 % data is a function returning the data that should be applied at the boundary.
202 % neighbour_scheme is an instance of Scheme that should be interfaced to. 210 % neighbour_scheme is an instance of Scheme that should be interfaced to.
203 % neighbour_boundary is a string specifying which boundary to interface to. 211 % neighbour_boundary is a string specifying which boundary to interface to.
204 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 212 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
205 default_arg('type','free'); 213 default_arg('type','free');
206 default_arg('parameter', []); 214 default_arg('parameter', []);
207
208 delta = @kroneckerDelta; % Kronecker delta
209 delta_b = @(i,j) 1-delta(i,j); % Logical not of Kronecker delta
210 215
211 % j is the coordinate direction of the boundary 216 % j is the coordinate direction of the boundary
212 % nj: outward unit normal component. 217 % nj: outward unit normal component.
213 % nj = -1 for west, south, bottom boundaries 218 % nj = -1 for west, south, bottom boundaries
214 % nj = 1 for east, north, top boundaries 219 % nj = 1 for east, north, top boundaries
224 229
225 E = obj.E; 230 E = obj.E;
226 Hi = obj.Hi; 231 Hi = obj.Hi;
227 H_gamma = obj.H_boundary{j}; 232 H_gamma = obj.H_boundary{j};
228 A = obj.A; 233 A = obj.A;
229 RHO = obj.RHO; 234 RHOi = obj.RHOi;
230 D1 = obj.D1;
231 235
232 phi = obj.phi; 236 phi = obj.phi;
233 H11 = obj.H11; 237 H11 = obj.H11;
234 h = obj.h; 238 h = obj.h;
239
240 RHOi_kron = obj.RHOi_kron;
241 Hi_kron = obj.Hi_kron;
242
243 % Divergence operator
244 Div = obj.Div{j};
235 245
236 switch type 246 switch type
237 % Dirichlet boundary condition 247 % Dirichlet boundary condition
238 case {'D','d','dirichlet'} 248 case {'D','d','dirichlet'}
239 tuning = 1.2; 249 tuning = 1.2;
240 phi = obj.phi; 250 phi = obj.phi;
241 251
242 % Initialize with zeros 252 sigma = tuning * obj.dim/(H11*h(j)) +...
243 m_tot = obj.grid.N(); 253 tuning * 1/(H11*h(j)*phi);
244 closure = sparse(m_tot*obj.dim, m_tot*obj.dim); 254
245 penalty = 0*E{j}*e{j}; 255 closure = - sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma*e{j}'*E{j}' ...
246 256 + nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma*e{j}'*E{j}';
247 % Loop over components to put penalties on 257
248 for i = 1:obj.dim 258 penalty = + sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma ...
249 sigma_i = tuning * obj.dim/(H11*h(j)) +... 259 - nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma;
250 tuning * 1/(H11*h(j)*phi);
251
252 closure = closure + E{i}*inv(RHO)*nj*Hi*...
253 ( ...
254 delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ...
255 + delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{j}' ...
256 ) ...
257 - delta(i,j)*sigma_i*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{j}';
258
259 penalty = penalty - ...
260 E{i}*inv(RHO)*nj*Hi*...
261 ( ...
262 delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ...
263 + delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' ...
264 ) ...
265 + delta(i,j)*sigma_i*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma;
266
267 end
268 260
269 % Free boundary condition 261 % Free boundary condition
270 case {'F','f','Free','free'} 262 case {'F','f','Free','free'}
271 closures = cell(obj.dim,1); 263
272 penalties = cell(obj.dim,obj.dim); 264 closure = -nj*E{j}*RHOi*Hi*e{j} ...
273
274 % Divergence operator
275 Div = obj.Div{j};
276
277 closure = -nj*E{j}*inv(RHO)*Hi*e{j} ...
278 *H_gamma*e{j}'*A*e{j}*e{j}'*Div; 265 *H_gamma*e{j}'*A*e{j}*e{j}'*Div;
279 penalty = nj*E{j}*inv(RHO)*Hi*e{j} ... 266 penalty = nj*E{j}*RHOi*Hi*e{j} ...
280 *H_gamma*e{j}'*A*e{j}; 267 *H_gamma*e{j}'*A*e{j};
281 268
282 269
283 % Unknown boundary condition 270 % Unknown boundary condition
284 otherwise 271 otherwise