Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 1059:d8ab528f10f6 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Massive boundary operator cleanup in Elastic2dVariable
author Martin Almquist <malmquist@stanford.edu>
date Sat, 26 Jan 2019 15:38:58 -0800
parents 84933722ec0e
children e40899094f20
comparison
equal deleted inserted replaced
1058:84933722ec0e 1059:d8ab528f10f6
13 dim 13 dim
14 14
15 order % Order of accuracy for the approximation 15 order % Order of accuracy for the approximation
16 16
17 % Diagonal matrices for varible coefficients 17 % Diagonal matrices for varible coefficients
18 LAMBDA % Variable coefficient, related to dilation 18 LAMBDA % Lame's first parameter, related to dilation
19 MU % Shear modulus, variable coefficient 19 MU % Shear modulus
20 RHO, RHOi % Density, variable 20 RHO, RHOi, RHOi_kron % Density
21 21
22 D % Total operator 22 D % Total operator
23 D1 % First derivatives 23 D1 % First derivatives
24 24
25 % Second derivatives 25 % Second derivatives
26 D2_lambda 26 D2_lambda
27 D2_mu 27 D2_mu
28 28
29 % Traction operators used for BC 29 % Boundary operators in cell format, used for BC
30 T_l, T_r 30 T_w, T_e, T_s, T_n
31 tau_l, tau_r 31
32 32 % Traction operators
33 H, Hi, H_1D % Inner products 33 tau_w, tau_e, tau_s, tau_n % Return vector field
34 e_l, e_r 34 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
35 35 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
36 36
37 d1_l, d1_r % Normal derivatives at the boundary 37 % Inner products
38 E % E{i}^T picks out component i 38 H, Hi, Hi_kron, H_1D
39 39
40 H_boundary % Boundary inner products 40 % Boundary inner products (for scalar field)
41 41 H_w, H_e, H_s, H_n
42 % Kroneckered norms and coefficients 42
43 RHOi_kron 43 % Boundary restriction operators
44 Hi_kron 44 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
45 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
46 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
47 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
48
49 % E{i}^T picks out component i
50 E
45 51
46 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. 52 % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
47 theta_R % Borrowing (d1- D1)^2 from R 53 theta_R % Borrowing (d1- D1)^2 from R
48 theta_H % First entry in norm matrix 54 theta_H % First entry in norm matrix
49 theta_M % Borrowing d1^2 from M. 55 theta_M % Borrowing d1^2 from M.
103 I = cell(dim,1); 109 I = cell(dim,1);
104 D1 = cell(dim,1); 110 D1 = cell(dim,1);
105 D2 = cell(dim,1); 111 D2 = cell(dim,1);
106 H = cell(dim,1); 112 H = cell(dim,1);
107 Hi = cell(dim,1); 113 Hi = cell(dim,1);
108 e_l = cell(dim,1); 114 e_0 = cell(dim,1);
109 e_r = cell(dim,1); 115 e_m = cell(dim,1);
110 d1_l = cell(dim,1); 116 d1_0 = cell(dim,1);
111 d1_r = cell(dim,1); 117 d1_m = cell(dim,1);
112 118
113 for i = 1:dim 119 for i = 1:dim
114 I{i} = speye(m(i)); 120 I{i} = speye(m(i));
115 D1{i} = ops{i}.D1; 121 D1{i} = ops{i}.D1;
116 D2{i} = ops{i}.D2; 122 D2{i} = ops{i}.D2;
117 H{i} = ops{i}.H; 123 H{i} = ops{i}.H;
118 Hi{i} = ops{i}.HI; 124 Hi{i} = ops{i}.HI;
119 e_l{i} = ops{i}.e_l; 125 e_0{i} = ops{i}.e_l;
120 e_r{i} = ops{i}.e_r; 126 e_m{i} = ops{i}.e_r;
121 d1_l{i} = ops{i}.d1_l; 127 d1_0{i} = ops{i}.d1_l;
122 d1_r{i} = ops{i}.d1_r; 128 d1_m{i} = ops{i}.d1_r;
123 end 129 end
124 130
125 %====== Assemble full operators ======== 131 %====== Assemble full operators ========
126 LAMBDA = spdiag(lambda); 132 LAMBDA = spdiag(lambda);
127 obj.LAMBDA = LAMBDA; 133 obj.LAMBDA = LAMBDA;
132 obj.RHOi = inv(RHO); 138 obj.RHOi = inv(RHO);
133 139
134 obj.D1 = cell(dim,1); 140 obj.D1 = cell(dim,1);
135 obj.D2_lambda = cell(dim,1); 141 obj.D2_lambda = cell(dim,1);
136 obj.D2_mu = cell(dim,1); 142 obj.D2_mu = cell(dim,1);
137 obj.e_l = cell(dim,1);
138 obj.e_r = cell(dim,1);
139 obj.d1_l = cell(dim,1);
140 obj.d1_r = cell(dim,1);
141 143
142 % D1 144 % D1
143 obj.D1{1} = kron(D1{1},I{2}); 145 obj.D1{1} = kron(D1{1},I{2});
144 obj.D1{2} = kron(I{1},D1{2}); 146 obj.D1{2} = kron(I{1},D1{2});
145 147
146 % Boundary operators 148 % Boundary restriction operators
147 obj.e_l{1} = kron(e_l{1},I{2}); 149 e_l = cell(dim,1);
148 obj.e_l{2} = kron(I{1},e_l{2}); 150 e_r = cell(dim,1);
149 obj.e_r{1} = kron(e_r{1},I{2}); 151 e_l{1} = kron(e_0{1}, I{2});
150 obj.e_r{2} = kron(I{1},e_r{2}); 152 e_l{2} = kron(I{1}, e_0{2});
151 153 e_r{1} = kron(e_m{1}, I{2});
152 obj.d1_l{1} = kron(d1_l{1},I{2}); 154 e_r{2} = kron(I{1}, e_m{2});
153 obj.d1_l{2} = kron(I{1},d1_l{2}); 155
154 obj.d1_r{1} = kron(d1_r{1},I{2}); 156 e_scalar_w = e_l{1};
155 obj.d1_r{2} = kron(I{1},d1_r{2}); 157 e_scalar_e = e_r{1};
158 e_scalar_s = e_l{2};
159 e_scalar_n = e_r{2};
160
161 I_dim = speye(dim, dim);
162 e_w = kron(e_scalar_w, I_dim);
163 e_e = kron(e_scalar_e, I_dim);
164 e_s = kron(e_scalar_s, I_dim);
165 e_n = kron(e_scalar_n, I_dim);
166
167 % Boundary derivatives
168 d1_l = cell(dim,1);
169 d1_r = cell(dim,1);
170 d1_l{1} = kron(d1_0{1}, I{2});
171 d1_l{2} = kron(I{1}, d1_0{2});
172 d1_r{1} = kron(d1_m{1}, I{2});
173 d1_r{2} = kron(I{1}, d1_m{2});
174
175
176 % E{i}^T picks out component i.
177 E = cell(dim,1);
178 I = speye(m_tot,m_tot);
179 for i = 1:dim
180 e = sparse(dim,1);
181 e(i) = 1;
182 E{i} = kron(I,e);
183 end
184 obj.E = E;
185
186 e1_w = (e_scalar_w'*E{1}')';
187 e1_e = (e_scalar_e'*E{1}')';
188 e1_s = (e_scalar_s'*E{1}')';
189 e1_n = (e_scalar_n'*E{1}')';
190
191 e2_w = (e_scalar_w'*E{2}')';
192 e2_e = (e_scalar_e'*E{2}')';
193 e2_s = (e_scalar_s'*E{2}')';
194 e2_n = (e_scalar_n'*E{2}')';
195
156 196
157 % D2 197 % D2
158 for i = 1:dim 198 for i = 1:dim
159 obj.D2_lambda{i} = sparse(m_tot); 199 obj.D2_lambda{i} = sparse(m_tot);
160 obj.D2_mu{i} = sparse(m_tot); 200 obj.D2_mu{i} = sparse(m_tot);
180 end 220 end
181 221
182 % Quadratures 222 % Quadratures
183 obj.H = kron(H{1},H{2}); 223 obj.H = kron(H{1},H{2});
184 obj.Hi = inv(obj.H); 224 obj.Hi = inv(obj.H);
185 obj.H_boundary = cell(dim,1); 225 obj.H_w = H{2};
186 obj.H_boundary{1} = H{2}; 226 obj.H_e = H{2};
187 obj.H_boundary{2} = H{1}; 227 obj.H_s = H{1};
228 obj.H_n = H{1};
188 obj.H_1D = {H{1}, H{2}}; 229 obj.H_1D = {H{1}, H{2}};
189
190 % E{i}^T picks out component i.
191 E = cell(dim,1);
192 I = speye(m_tot,m_tot);
193 for i = 1:dim
194 e = sparse(dim,1);
195 e(i) = 1;
196 E{i} = kron(I,e);
197 end
198 obj.E = E;
199 230
200 % Differentiation matrix D (without SAT) 231 % Differentiation matrix D (without SAT)
201 D2_lambda = obj.D2_lambda; 232 D2_lambda = obj.D2_lambda;
202 D2_mu = obj.D2_mu; 233 D2_mu = obj.D2_mu;
203 D1 = obj.D1; 234 D1 = obj.D1;
219 %=========================================%' 250 %=========================================%'
220 251
221 % Numerical traction operators for BC. 252 % Numerical traction operators for BC.
222 % Because d1 =/= e0^T*D1, the numerical tractions are different 253 % Because d1 =/= e0^T*D1, the numerical tractions are different
223 % at every boundary. 254 % at every boundary.
255 %
256 % Formula at boundary j: % tau^{j}_i = sum_k T^{j}_{ik} u_k
257 %
224 T_l = cell(dim,1); 258 T_l = cell(dim,1);
225 T_r = cell(dim,1); 259 T_r = cell(dim,1);
226 tau_l = cell(dim,1); 260 tau_l = cell(dim,1);
227 tau_r = cell(dim,1); 261 tau_r = cell(dim,1);
228 % tau^{j}_i = sum_k T^{j}_{ik} u_k 262
229
230 d1_l = obj.d1_l;
231 d1_r = obj.d1_r;
232 e_l = obj.e_l;
233 e_r = obj.e_r;
234 D1 = obj.D1; 263 D1 = obj.D1;
235 264
236 % Loop over boundaries 265 % Loop over boundaries
237 for j = 1:dim 266 for j = 1:dim
238 T_l{j} = cell(dim,dim); 267 T_l{j} = cell(dim,dim);
248 [~, n_l] = size(e_l{j}); 277 [~, n_l] = size(e_l{j});
249 [~, n_r] = size(e_r{j}); 278 [~, n_r] = size(e_r{j});
250 279
251 % Loop over components 280 % Loop over components
252 for i = 1:dim 281 for i = 1:dim
253 tau_l{j}{i} = sparse(n_l, dim*m_tot); 282 tau_l{j}{i} = sparse(dim*m_tot, n_l);
254 tau_r{j}{i} = sparse(n_r, dim*m_tot); 283 tau_r{j}{i} = sparse(dim*m_tot, n_r);
255 for k = 1:dim 284 for k = 1:dim
256 T_l{j}{i,k} = ... 285 T_l{j}{i,k} = ...
257 -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})... 286 (-d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})...
258 -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})... 287 -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})...
259 -d(i,k)*MU_l*d1_l{j}'; 288 -d(i,k)*MU_l*d1_l{j}')';
260 289
261 T_r{j}{i,k} = ... 290 T_r{j}{i,k} = ...
262 d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})... 291 (d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})...
263 +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})... 292 +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})...
264 +d(i,k)*MU_r*d1_r{j}'; 293 +d(i,k)*MU_r*d1_r{j}')';
265 294
266 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; 295 tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,k}'*E{k}')';
267 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; 296 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,k}'*E{k}')';
268 end 297 end
269 298
270 end 299 end
271 end 300 end
272 301
273 % Transpose T and tau to match boundary operator convention 302 % Traction tensors, T_ij
274 for i = 1:dim 303 obj.T_w = T_l{1};
275 for j = 1:dim 304 obj.T_e = T_r{1};
276 tau_l{i}{j} = transpose(tau_l{i}{j}); 305 obj.T_s = T_l{2};
277 tau_r{i}{j} = transpose(tau_r{i}{j}); 306 obj.T_n = T_r{2};
278 for k = 1:dim 307
279 T_l{i}{j,k} = transpose(T_l{i}{j,k}); 308 % Restriction operators
280 T_r{i}{j,k} = transpose(T_r{i}{j,k}); 309 obj.e_w = e_w;
281 end 310 obj.e_e = e_e;
282 end 311 obj.e_s = e_s;
283 end 312 obj.e_n = e_n;
284 313
285 obj.T_l = T_l; 314 obj.e1_w = e1_w;
286 obj.T_r = T_r; 315 obj.e1_e = e1_e;
287 obj.tau_l = tau_l; 316 obj.e1_s = e1_s;
288 obj.tau_r = tau_r; 317 obj.e1_n = e1_n;
318
319 obj.e2_w = e2_w;
320 obj.e2_e = e2_e;
321 obj.e2_s = e2_s;
322 obj.e2_n = e2_n;
323
324 obj.e_scalar_w = e_scalar_w;
325 obj.e_scalar_e = e_scalar_e;
326 obj.e_scalar_s = e_scalar_s;
327 obj.e_scalar_n = e_scalar_n;
328
329 % First component of traction
330 obj.tau1_w = tau_l{1}{1};
331 obj.tau1_e = tau_r{1}{1};
332 obj.tau1_s = tau_l{2}{1};
333 obj.tau1_n = tau_r{2}{1};
334
335 % Second component of traction
336 obj.tau2_w = tau_l{1}{2};
337 obj.tau2_e = tau_r{1}{2};
338 obj.tau2_s = tau_l{2}{2};
339 obj.tau2_n = tau_r{2}{2};
340
341 % Traction vectors
342 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
343 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
344 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
345 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
289 346
290 % Kroneckered norms and coefficients 347 % Kroneckered norms and coefficients
291 I_dim = speye(dim);
292 obj.RHOi_kron = kron(obj.RHOi, I_dim); 348 obj.RHOi_kron = kron(obj.RHOi, I_dim);
293 obj.Hi_kron = kron(obj.Hi, I_dim); 349 obj.Hi_kron = kron(obj.Hi, I_dim);
294 350
295 % Misc. 351 % Misc.
296 obj.m = m; 352 obj.m = m;
358 type = bc{2}; 414 type = bc{2};
359 if ischar(comp) 415 if ischar(comp)
360 comp = obj.getComponent(comp, boundary); 416 comp = obj.getComponent(comp, boundary);
361 end 417 end
362 418
363 e = obj.getBoundaryOperator('e', boundary); 419 e = obj.getBoundaryOperatorForScalarField('e', boundary);
364 T = obj.getBoundaryOperator('T', boundary); 420 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
365 tau = obj.getBoundaryOperator('tau', boundary); 421 T = obj.getBoundaryTractionOperator(boundary);
366 H_gamma = obj.getBoundaryOperator('H', boundary); 422 alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary);
423 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
367 424
368 E = obj.E; 425 E = obj.E;
369 Hi = obj.Hi; 426 Hi = obj.Hi;
370 LAMBDA = obj.LAMBDA; 427 LAMBDA = obj.LAMBDA;
371 MU = obj.MU; 428 MU = obj.MU;
373 430
374 dim = obj.dim; 431 dim = obj.dim;
375 m_tot = obj.grid.N(); 432 m_tot = obj.grid.N();
376 433
377 % Preallocate 434 % Preallocate
378 [~, col] = size(tau{comp}); 435 [~, col] = size(tau);
379 closure = sparse(dim*m_tot, dim*m_tot); 436 closure = sparse(dim*m_tot, dim*m_tot);
380 penalty = sparse(dim*m_tot, col); 437 penalty = sparse(dim*m_tot, col);
381 438
382 k = comp; 439 k = comp;
383 switch type 440 switch type
384 441
385 % Dirichlet boundary condition 442 % Dirichlet boundary condition
386 case {'D','d','dirichlet','Dirichlet'} 443 case {'D','d','dirichlet','Dirichlet'}
387
388 alpha = obj.getBoundaryOperator('alpha', boundary);
389 444
390 % Loop over components that Dirichlet penalties end up on 445 % Loop over components that Dirichlet penalties end up on
391 for i = 1:dim 446 for i = 1:dim
392 C = transpose(T{k,i}); 447 C = transpose(T{k,i});
393 A = -tuning*e*transpose(alpha{i,k}); 448 A = -tuning*e*transpose(alpha{i,k});
396 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 451 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
397 end 452 end
398 453
399 % Free boundary condition 454 % Free boundary condition
400 case {'F','f','Free','free','traction','Traction','t','T'} 455 case {'F','f','Free','free','traction','Traction','t','T'}
401 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}'; 456 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau';
402 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; 457 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
403 458
404 % Unknown boundary condition 459 % Unknown boundary condition
405 otherwise 460 otherwise
406 error('No such boundary condition: type = %s',type); 461 error('No such boundary condition: type = %s',type);
433 % u denotes the solution in the own domain 488 % u denotes the solution in the own domain
434 % v denotes the solution in the neighbour domain 489 % v denotes the solution in the neighbour domain
435 % Operators without subscripts are from the own domain. 490 % Operators without subscripts are from the own domain.
436 491
437 % Get boundary operators 492 % Get boundary operators
438 e = obj.getBoundaryOperator('e_tot', boundary); 493 e = obj.getBoundaryOperator('e', boundary);
439 tau = obj.getBoundaryOperator('tau_tot', boundary); 494 tau = obj.getBoundaryOperator('tau', boundary);
440 495
441 e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary); 496 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
442 tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary); 497 tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary);
443 498
444 H_gamma = obj.getBoundaryQuadrature(boundary); 499 H_gamma = obj.getBoundaryQuadrature(boundary);
445 500
446 % Operators and quantities that correspond to the own domain only 501 % Operators and quantities that correspond to the own domain only
447 Hi = obj.Hi_kron; 502 Hi = obj.Hi_kron;
448 RHOi = obj.RHOi_kron; 503 RHOi = obj.RHOi_kron;
449 504
450 % Penalty strength operators 505 % Penalty strength operators
451 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary); 506 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary);
452 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary); 507 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary);
453 508
454 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e'); 509 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
455 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v'); 510 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
456 511
457 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau'; 512 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
490 end 545 end
491 end 546 end
492 547
493 % Returns the boundary operator op for the boundary specified by the string boundary. 548 % Returns the boundary operator op for the boundary specified by the string boundary.
494 % op -- string 549 % op -- string
495 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
496 function o = getBoundaryOperator(obj, op, boundary) 550 function o = getBoundaryOperator(obj, op, boundary)
497 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 551 assertIsMember(boundary, {'w', 'e', 's', 'n'})
498 assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'}) 552 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha'})
499
500 switch boundary
501 case {'w', 'e'}
502 j = 1;
503 case {'s', 'n'}
504 j = 2;
505 end
506 553
507 switch op 554 switch op
508 case 'e' 555
509 switch boundary 556 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
510 case {'w', 's'} 557 o = obj.([op, '_', boundary]);
511 o = obj.e_l{j};
512 case {'e', 'n'}
513 o = obj.e_r{j};
514 end
515
516 case 'e_tot'
517 e = obj.getBoundaryOperator('e', boundary);
518 I_dim = speye(obj.dim, obj.dim);
519 o = kron(e, I_dim);
520
521 case 'd'
522 switch boundary
523 case {'w', 's'}
524 o = obj.d1_l{j};
525 case {'e', 'n'}
526 o = obj.d1_r{j};
527 end
528
529 case 'T'
530 switch boundary
531 case {'w', 's'}
532 o = obj.T_l{j};
533 case {'e', 'n'}
534 o = obj.T_r{j};
535 end
536
537 case 'tau'
538 switch boundary
539 case {'w', 's'}
540 o = obj.tau_l{j};
541 case {'e', 'n'}
542 o = obj.tau_r{j};
543 end
544
545 case 'tau_tot'
546 e = obj.getBoundaryOperator('e', boundary);
547 tau = obj.getBoundaryOperator('tau', boundary);
548
549 I_dim = speye(obj.dim, obj.dim);
550 e_tot = kron(e, I_dim);
551 E = obj.E;
552 tau_tot = (e_tot'*E{1}*e*tau{1}')';
553 for i = 2:obj.dim
554 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
555 end
556 o = tau_tot;
557
558 case 'H'
559 o = obj.H_boundary{j};
560 558
561 case 'alpha' 559 case 'alpha'
562 % alpha = alpha(i,j) is the penalty strength for displacement BC. 560 % alpha = alpha(i,j) is the penalty strength for displacement BC.
563 e = obj.getBoundaryOperator('e', boundary); 561 e = obj.getBoundaryOperatorForScalarField('e', boundary);
564 562 e_tot = obj.getBoundaryOperator('e', boundary);
565 LAMBDA = obj.LAMBDA; 563 alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary);
566 MU = obj.MU;
567
568 dim = obj.dim;
569 theta_R = obj.theta_R{j};
570 theta_H = obj.theta_H{j};
571 theta_M = obj.theta_M{j};
572
573 a_lambda = dim/theta_H + 1/theta_R;
574 a_mu_i = 2/theta_M;
575 a_mu_ij = 2/theta_H + 1/theta_R;
576
577 d = @kroneckerDelta; % Kronecker delta
578 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
579 alpha = cell(obj.dim, obj.dim);
580
581 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
582 + d(i,j)* a_mu_i*MU ...
583 + db(i,j)*a_mu_ij*MU;
584 for i = 1:obj.dim
585 for l = 1:obj.dim
586 alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
587 end
588 end
589
590 o = alpha;
591
592 case 'alpha_tot'
593 % alpha = alpha(i,j) is the penalty strength for displacement BC.
594 e = obj.getBoundaryOperator('e', boundary);
595 e_tot = obj.getBoundaryOperator('e_tot', boundary);
596 alpha = obj.getBoundaryOperator('alpha', boundary);
597 E = obj.E; 564 E = obj.E;
598 [m, n] = size(alpha{1,1}); 565 [m, n] = size(alpha{1,1});
599 alpha_tot = sparse(m*obj.dim, n*obj.dim); 566 alpha_tot = sparse(m*obj.dim, n*obj.dim);
600 for i = 1:obj.dim 567 for i = 1:obj.dim
601 for l = 1:obj.dim 568 for l = 1:obj.dim
605 o = alpha_tot; 572 o = alpha_tot;
606 end 573 end
607 574
608 end 575 end
609 576
577 % Returns the boundary operator op for the boundary specified by the string boundary.
578 % op -- string
579 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
580 assertIsMember(boundary, {'w', 'e', 's', 'n'})
581 assertIsMember(op, {'e', 'alpha'})
582
583 switch op
584
585 case 'e'
586 o = obj.(['e_scalar', '_', boundary]);
587
588 case 'alpha'
589 % alpha = alpha(i,j) is the penalty strength for displacement BC.
590 e = obj.getBoundaryOperatorForScalarField('e', boundary);
591
592 LAMBDA = obj.LAMBDA;
593 MU = obj.MU;
594 dim = obj.dim;
595
596 switch boundary
597 case {'w', 'e'}
598 k = 1;
599 case {'s', 'n'}
600 k = 2;
601 end
602
603 theta_R = obj.theta_R{k};
604 theta_H = obj.theta_H{k};
605 theta_M = obj.theta_M{k};
606
607 a_lambda = dim/theta_H + 1/theta_R;
608 a_mu_i = 2/theta_M;
609 a_mu_ij = 2/theta_H + 1/theta_R;
610
611 d = @kroneckerDelta; % Kronecker delta
612 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
613 alpha = cell(obj.dim, obj.dim);
614
615 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
616 + d(i,j)* a_mu_i*MU ...
617 + db(i,j)*a_mu_ij*MU;
618 for i = 1:obj.dim
619 for j = 1:obj.dim
620 alpha{i,j} = d(i,j)*alpha_func(i,k)*e;
621 end
622 end
623
624 o = alpha;
625 end
626
627 end
628
629 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
630 % op -- string
631 function T = getBoundaryTractionOperator(obj, boundary)
632 assertIsMember(boundary, {'w', 'e', 's', 'n'})
633
634 T = obj.(['T', '_', boundary]);
635 end
636
637 % Returns square boundary quadrature matrix, of dimension
638 % corresponding to the number of boundary unknowns
639 %
640 % boundary -- string
610 function H = getBoundaryQuadrature(obj, boundary) 641 function H = getBoundaryQuadrature(obj, boundary)
611 switch boundary 642 assertIsMember(boundary, {'w', 'e', 's', 'n'})
612 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 643
613 j = 1; 644 H = obj.getBoundaryQuadratureForScalarField(boundary);
614 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
615 j = 2;
616 otherwise
617 error('No such boundary: boundary = %s',boundary);
618 end
619 H = obj.H_boundary{j};
620 I_dim = speye(obj.dim, obj.dim); 645 I_dim = speye(obj.dim, obj.dim);
621 H = kron(H, I_dim); 646 H = kron(H, I_dim);
622 end 647 end
623 648
649 % Returns square boundary quadrature matrix, of dimension
650 % corresponding to the number of boundary grid points
651 %
652 % boundary -- string
653 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
654 assertIsMember(boundary, {'w', 'e', 's', 'n'})
655
656 H_b = obj.(['H_', boundary]);
657 end
658
624 function N = size(obj) 659 function N = size(obj)
625 N = obj.dim*prod(obj.m); 660 N = obj.dim*prod(obj.m);
626 end 661 end
627 end 662 end
628 end 663 end