Mercurial > repos > public > sbplib

comparison +scheme/elasticVariable.m @ 678:06676c40e77f feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add scheme for complete elastic operator. Traction BC working with MMS variable coefficient.
author Martin Almquist <malmquist@stanford.edu>
date Mon, 05 Feb 2018 11:06:15 -0800
parents
children 247b58a4dbe8
comparison
equal deleted inserted replaced
677:eeaf9a00e304 678:06676c40e77f
1 classdef elasticVariable < scheme.Scheme
2
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
5
6
7 properties
8 m % Number of points in each direction, possibly a vector
9 h % Grid spacing
10
11 grid
12 dim
13
14 order % Order of accuracy for the approximation
15
16 % Diagonal matrices for varible coefficients
17 LAMBDA % Variable coefficient, related to dilation
18 MU % Shear modulus, variable coefficient
19 RHO, RHOi % Density
20
21 D % Total operator
22 D1 % First derivatives
23
24 % Second derivatives
25 D2_lambda
26 D2_mu
27
28 % Traction operators used for BC
29 T_l, T_r
30 tau_l, tau_r
31
32 H, Hi % Inner products
33 phi % Borrowing constant for (d1 - e^T*D1) from R
34 H11 % First element of H
35 e_l, e_r
36 d1_l, d1_r % Normal derivatives at the boundary
37 E % E{i}^T picks out component i
38
39 H_boundary % Boundary inner products
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
47 RHOi_kron
48 Hi_kron
49 end
50
51 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
56 function obj = elasticVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
57 default_arg('opSet',@sbp.D2Variable);
58 default_arg('lambda_fun', @(x,y) 0*x+1);
59 default_arg('mu_fun', @(x,y) 0*x+1);
60 default_arg('rho_fun', @(x,y) 0*x+1);
61 dim = 2;
62
63 assert(isa(g, 'grid.Cartesian'))
64
65 lambda = grid.evalOn(g, lambda_fun);
66 mu = grid.evalOn(g, mu_fun);
67 rho = grid.evalOn(g, rho_fun);
68 m = g.size();
69 m_tot = g.N();
70
71 h = g.scaling();
72 L = (m-1).*h;
73
74 % 1D operators
75 ops = cell(dim,1);
76 for i = 1:dim
77 ops{i} = opSet(m(i), {0, L(i)}, order);
78 end
79
80 % Borrowing constants
81 beta = ops{1}.borrowing.R.delta_D;
82 obj.H11 = ops{1}.borrowing.H11;
83 obj.phi = beta/obj.H11;
84
85 I = cell(dim,1);
86 D1 = cell(dim,1);
87 D2 = cell(dim,1);
88 H = cell(dim,1);
89 Hi = cell(dim,1);
90 e_l = cell(dim,1);
91 e_r = cell(dim,1);
92 d1_l = cell(dim,1);
93 d1_r = cell(dim,1);
94
95 for i = 1:dim
96 I{i} = speye(m(i));
97 D1{i} = ops{i}.D1;
98 D2{i} = ops{i}.D2;
99 H{i} = ops{i}.H;
100 Hi{i} = ops{i}.HI;
101 e_l{i} = ops{i}.e_l;
102 e_r{i} = ops{i}.e_r;
103 d1_l{i} = ops{i}.d1_l;
104 d1_r{i} = ops{i}.d1_r;
105 end
106
107 %====== Assemble full operators ========
108 LAMBDA = spdiag(lambda);
109 obj.LAMBDA = LAMBDA;
110 MU = spdiag(mu);
111 obj.MU = MU;
112 RHO = spdiag(rho);
113 obj.RHO = RHO;
114 obj.RHOi = inv(RHO);
115
116 obj.D1 = cell(dim,1);
117 obj.D2_lambda = cell(dim,1);
118 obj.D2_mu = cell(dim,1);
119 obj.e_l = cell(dim,1);
120 obj.e_r = cell(dim,1);
121 obj.d1_l = cell(dim,1);
122 obj.d1_r = cell(dim,1);
123
124 % D1
125 obj.D1{1} = kron(D1{1},I{2});
126 obj.D1{2} = kron(I{1},D1{2});
127
128 % Boundary operators
129 obj.e_l{1} = kron(e_l{1},I{2});
130 obj.e_l{2} = kron(I{1},e_l{2});
131 obj.e_r{1} = kron(e_r{1},I{2});
132 obj.e_r{2} = kron(I{1},e_r{2});
133
134 obj.d1_l{1} = kron(d1_l{1},I{2});
135 obj.d1_l{2} = kron(I{1},d1_l{2});
136 obj.d1_r{1} = kron(d1_r{1},I{2});
137 obj.d1_r{2} = kron(I{1},d1_r{2});
138
139 % D2
140 for i = 1:dim
141 obj.D2_lambda{i} = sparse(m_tot);
142 obj.D2_mu{i} = sparse(m_tot);
143 end
144 ind = grid.funcToMatrix(g, 1:m_tot);
145
146 for i = 1:m(2)
147 D_lambda = D2{1}(lambda(ind(:,i)));
148 D_mu = D2{1}(mu(ind(:,i)));
149
150 p = ind(:,i);
151 obj.D2_lambda{1}(p,p) = D_lambda;
152 obj.D2_mu{1}(p,p) = D_mu;
153 end
154
155 for i = 1:m(1)
156 D_lambda = D2{2}(lambda(ind(i,:)));
157 D_mu = D2{2}(mu(ind(i,:)));
158
159 p = ind(i,:);
160 obj.D2_lambda{2}(p,p) = D_lambda;
161 obj.D2_mu{2}(p,p) = D_mu;
162 end
163
164 % Quadratures
165 obj.H = kron(H{1},H{2});
166 obj.Hi = inv(obj.H);
167 obj.H_boundary = cell(dim,1);
168 obj.H_boundary{1} = H{2};
169 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
183 % E{i}^T picks out component i.
184 E = cell(dim,1);
185 I = speye(m_tot,m_tot);
186 for i = 1:dim
187 e = sparse(dim,1);
188 e(i) = 1;
189 E{i} = kron(I,e);
190 end
191 obj.E = E;
192
193 % Differentiation matrix D (without SAT)
194 D2_lambda = obj.D2_lambda;
195 D2_mu = obj.D2_mu;
196 D1 = obj.D1;
197 D = sparse(dim*m_tot,dim*m_tot);
198 d = @kroneckerDelta; % Kronecker delta
199 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
200 for i = 1:dim
201 for j = 1:dim
202 D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +...
203 db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ...
204 );
205 D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +...
206 db(i,j)*D1{j}*MU*D1{i}*E{j}' + ...
207 D2_mu{j}*E{i}' ...
208 );
209 end
210 end
211 obj.D = D;
212 %=========================================%
213
214 % Numerical traction operators for BC.
215 % Because d1 =/= e0^T*D1, the numerical tractions are different
216 % at every boundary.
217 T_l = cell(dim,1);
218 T_r = cell(dim,1);
219 tau_l = cell(dim,1);
220 tau_r = cell(dim,1);
221 % tau^{j}_i = sum_k T^{j}_{ik} u_k
222
223 d1_l = obj.d1_l;
224 d1_r = obj.d1_r;
225 e_l = obj.e_l;
226 e_r = obj.e_r;
227 D1 = obj.D1;
228
229 % Loop over boundaries
230 for j = 1:dim
231 T_l{j} = cell(dim,dim);
232 T_r{j} = cell(dim,dim);
233 tau_l{j} = cell(dim,1);
234 tau_r{j} = cell(dim,1);
235
236 % Loop over components
237 for i = 1:dim
238 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
239 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
240 for k = 1:dim
241 T_l{j}{i,k} = ...
242 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
243 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
244 -d(i,k)*MU*e_l{j}*d1_l{j}';
245
246 T_r{j}{i,k} = ...
247 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
248 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
249 +d(i,k)*MU*e_r{j}*d1_r{j}';
250
251 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
252 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
253 end
254
255 end
256 end
257 obj.T_l = T_l;
258 obj.T_r = T_r;
259 obj.tau_l = tau_l;
260 obj.tau_r = tau_r;
261
262 % Kroneckered norms and coefficients
263 I_dim = speye(dim);
264 obj.RHOi_kron = kron(obj.RHOi, I_dim);
265 obj.Hi_kron = kron(obj.Hi, I_dim);
266
267 % Misc.
268 obj.m = m;
269 obj.h = h;
270 obj.order = order;
271 obj.grid = g;
272 obj.dim = dim;
273
274 end
275
276
277 % 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.
279 % 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.
281 % 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.
283 % neighbour_boundary is a string specifying which boundary to interface to.
284 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
285 default_arg('type','free');
286 default_arg('parameter', []);
287
288 % j is the coordinate direction of the boundary
289 % nj: outward unit normal component.
290 % nj = -1 for west, south, bottom boundaries
291 % nj = 1 for east, north, top boundaries
292 [j, nj] = obj.get_boundary_number(boundary);
293 switch nj
294 case 1
295 e = obj.e_r;
296 d = obj.d1_r;
297 tau = obj.tau_r{j};
298 T = obj.T_r{j};
299 case -1
300 e = obj.e_l;
301 d = obj.d1_l;
302 tau = obj.tau_l{j};
303 T = obj.T_l{j};
304 end
305
306 E = obj.E;
307 Hi = obj.Hi;
308 H_gamma = obj.H_boundary{j};
309 LAMBDA = obj.LAMBDA;
310 MU = obj.MU;
311 RHOi = obj.RHOi;
312
313 phi = obj.phi;
314 H11 = obj.H11;
315 h = obj.h;
316 dim = obj.dim;
317 m_tot = obj.grid.N();
318
319 RHOi_kron = obj.RHOi_kron;
320 Hi_kron = obj.Hi_kron;
321
322 closure = sparse(dim*m_tot, dim*m_tot);
323 penalty = cell(dim,1);
324 switch type
325 % Dirichlet boundary condition
326 case {'D','d','dirichlet'}
327 error('not implemented')
328 tuning = 1.2;
329 phi = obj.phi;
330
331 sigma = tuning * obj.dim/(H11*h(j)) +...
332 tuning * 1/(H11*h(j)*phi);
333
334 closure = - sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma*e{j}'*E{j}' ...
335 + nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma*e{j}'*E{j}';
336
337 penalty = + sigma*E{j}*RHOi*Hi*A*e{j}*H_gamma ...
338 - nj*RHOi_kron*Hi_kron*Div'*A*e{j}*H_gamma;
339
340 % Free boundary condition
341 case {'F','f','Free','free'}
342
343 % Loop over components of traction
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
350 % Unknown boundary condition
351 otherwise
352 error('No such boundary condition: type = %s',type);
353 end
354 end
355
356 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
357 % u denotes the solution in the own domain
358 % v denotes the solution in the neighbour domain
359 tuning = 1.2;
360 % tuning = 20.2;
361 error('Interface not implemented');
362 end
363
364 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
365 function [j, nj] = get_boundary_number(obj, boundary)
366
367 switch boundary
368 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
369 j = 1;
370 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
371 j = 2;
372 otherwise
373 error('No such boundary: boundary = %s',boundary);
374 end
375
376 switch boundary
377 case {'w','W','west','West','s','S','south','South'}
378 nj = -1;
379 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
380 nj = 1;
381 end
382 end
383
384 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
385 function [return_op] = get_boundary_operator(obj, op, boundary)
386
387 switch boundary
388 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
389 j = 1;
390 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
391 j = 2;
392 otherwise
393 error('No such boundary: boundary = %s',boundary);
394 end
395
396 switch op
397 case 'e'
398 switch boundary
399 case {'w','W','west','West','s','S','south','South'}
400 return_op = obj.e_l{j};
401 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
402 return_op = obj.e_r{j};
403 end
404 case 'd'
405 switch boundary
406 case {'w','W','west','West','s','S','south','South'}
407 return_op = obj.d_l{j};
408 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
409 return_op = obj.d_r{j};
410 end
411 otherwise
412 error(['No such operator: operatr = ' op]);
413 end
414
415 end
416
417 function N = size(obj)
418 N = prod(obj.m);
419 end
420 end
421 end