Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dStaggeredCurvilinearAnisotropic.m @ 1281:01a0500de446 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add working implementation of Elastic2dStaggeredCurvilinearAnisotropic, which wraps around the Cartesian anisotropic scheme.
author Martin Almquist <malmquist@stanford.edu>
date Sun, 07 Jun 2020 11:33:44 -0700
parents
children
comparison
equal deleted inserted replaced
1280:c7f5b480e455 1281:01a0500de446
1 classdef Elastic2dStaggeredCurvilinearAnisotropic < scheme.Scheme
2
3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
5 % in curvilinear coordinates, using Lebedev staggered grids
6
7 properties
8 m % Number of points in each direction, possibly a vector
9 h % Grid spacing
10
11 grid
12 dim
13 nGrids
14
15 order % Order of accuracy for the approximation
16
17 % Diagonal matrices for variable coefficients
18 J, Ji
19 RHO % Density
20 C % Elastic stiffness tensor
21
22 D % Total operator
23
24 % Dx, Dy % Physical derivatives
25 n_w, n_e, n_s, n_n % Physical normals
26
27 % Boundary operators in cell format, used for BC
28 T_w, T_e, T_s, T_n
29
30 % Traction operators
31 % tau_w, tau_e, tau_s, tau_n % Return vector field
32 % tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
33 % tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
34
35 % Inner products
36 H
37
38 % Boundary inner products (for scalar field)
39 H_w, H_e, H_s, H_n
40
41 % Surface Jacobian vectors
42 s_w, s_e, s_s, s_n
43
44 % Boundary restriction operators
45 e_w_u, e_e_u, e_s_u, e_n_u % Act on scalar field, return scalar field at boundary
46 e_w_s, e_e_s, e_s_s, e_n_s % Act on scalar field, return scalar field at boundary
47 % e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
48 % e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
49 % e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
50 % en_w, en_e, en_s, en_n % Act on vector field, return normal component
51
52 % U{i}^T picks out component i
53 U
54
55 % G{i}^T picks out displacement grid i
56 G
57
58 % Elastic2dVariableAnisotropic object for reference domain
59 refObj
60 end
61
62 methods
63
64 % The coefficients can either be function handles or grid functions
65 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
66 function obj = Elastic2dStaggeredCurvilinearAnisotropic(g, order, rho, C)
67 default_arg('rho', @(x,y) 0*x+1);
68
69 opSet = @sbp.D1StaggeredUpwind;
70 dim = 2;
71 nGrids = 2;
72
73 C_default = cell(dim,dim,dim,dim);
74 for i = 1:dim
75 for j = 1:dim
76 for k = 1:dim
77 for l = 1:dim
78 C_default{i,j,k,l} = @(x,y) 0*x;
79 end
80 end
81 end
82 end
83 default_arg('C', C_default);
84 assert(isa(g, 'grid.Staggered'));
85
86 g_u = g.gridGroups{1};
87 g_s = g.gridGroups{2};
88
89 m_u = {g_u{1}.size(), g_u{2}.size()};
90 m_s = {g_s{1}.size(), g_s{2}.size()};
91
92 if isa(rho, 'function_handle')
93 rho_vec = cell(nGrids, 1);
94 for a = 1:nGrids
95 rho_vec{a} = grid.evalOn(g_u{a}, rho);
96 end
97 rho = rho_vec;
98 end
99 for a = 1:nGrids
100 RHO{a} = spdiag(rho{a});
101 end
102 obj.RHO = RHO;
103
104 C_mat = cell(nGrids, 1);
105 for a = 1:nGrids
106 C_mat{a} = cell(dim,dim,dim,dim);
107 end
108 for a = 1:nGrids
109 for i = 1:dim
110 for j = 1:dim
111 for k = 1:dim
112 for l = 1:dim
113 if numel(C) == dim
114 C_mat{a}{i,j,k,l} = spdiag(C{a}{i,j,k,l});
115 else
116 C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
117 end
118 end
119 end
120 end
121 end
122 end
123 obj.C = C_mat;
124
125 C = cell(nGrids, 1);
126 for a = 1:nGrids
127 C{a} = cell(dim,dim,dim,dim);
128 for i = 1:dim
129 for j = 1:dim
130 for k = 1:dim
131 for l = 1:dim
132 C{a}{i,j,k,l} = diag(C_mat{a}{i,j,k,l});
133 end
134 end
135 end
136 end
137 end
138
139 % Reference m for primal grid
140 m = g_u{1}.size();
141
142 % 1D operators
143 ops = cell(dim,1);
144 D1p = cell(dim, 1);
145 D1d = cell(dim, 1);
146 mp = cell(dim, 1);
147 md = cell(dim, 1);
148 Ip = cell(dim, 1);
149 Id = cell(dim, 1);
150 Hp = cell(dim, 1);
151 Hd = cell(dim, 1);
152
153 opSet = @sbp.D1StaggeredUpwind;
154 for i = 1:dim
155 ops{i} = opSet(m(i), {0,1}, order);
156 D1p{i} = ops{i}.D1_dual;
157 D1d{i} = ops{i}.D1_primal;
158 mp{i} = length(ops{i}.x_primal);
159 md{i} = length(ops{i}.x_dual);
160 Ip{i} = speye(mp{i}, mp{i});
161 Id{i} = speye(md{i}, md{i});
162 Hp{i} = ops{i}.H_primal;
163 Hd{i} = ops{i}.H_dual;
164 ep_l{i} = ops{i}.e_primal_l;
165 ep_r{i} = ops{i}.e_primal_r;
166 ed_l{i} = ops{i}.e_dual_l;
167 ed_r{i} = ops{i}.e_dual_r;
168 end
169
170 % D1_u2s{a, b}{i} approximates ddi and
171 % takes from u grid number b to s grid number a
172 % Some of D1_x2y{a, b} are 0.
173 D1_u2s = cell(nGrids, nGrids);
174 D1_s2u = cell(nGrids, nGrids);
175
176 N_u = cell(nGrids, 1);
177 N_s = cell(nGrids, 1);
178 for a = 1:nGrids
179 N_u{a} = g_u{a}.N();
180 N_s{a} = g_s{a}.N();
181 end
182
183 %---- Grid layout -------
184 % gu1 = xp o yp;
185 % gu2 = xd o yd;
186 % gs1 = xd o yp;
187 % gs2 = xp o yd;
188 %------------------------
189
190 % Logical operators
191 D1_u2s{1,1}{1} = kron(D1p{1}, Ip{2});
192 D1_s2u{1,1}{1} = kron(D1d{1}, Ip{2});
193
194 D1_u2s{1,2}{2} = kron(Id{1}, D1d{2});
195 D1_u2s{2,1}{2} = kron(Ip{1}, D1p{2});
196
197 D1_s2u{1,2}{2} = kron(Ip{1}, D1d{2});
198 D1_s2u{2,1}{2} = kron(Id{1}, D1p{2});
199
200 D1_u2s{2,2}{1} = kron(D1d{1}, Id{2});
201 D1_s2u{2,2}{1} = kron(D1p{1}, Id{2});
202
203 D1_u2s{1,1}{2} = sparse(N_s{1}, N_u{1});
204 D1_s2u{1,1}{2} = sparse(N_u{1}, N_s{1});
205
206 D1_u2s{2,2}{2} = sparse(N_s{2}, N_u{2});
207 D1_s2u{2,2}{2} = sparse(N_u{2}, N_s{2});
208
209 D1_u2s{1,2}{1} = sparse(N_s{1}, N_u{2});
210 D1_s2u{1,2}{1} = sparse(N_u{1}, N_s{2});
211
212 D1_u2s{2,1}{1} = sparse(N_s{2}, N_u{1});
213 D1_s2u{2,1}{1} = sparse(N_u{2}, N_s{1});
214
215
216 %---- Combine grids and components -----
217
218 % U{a}{i}^T picks out u component i on grid a
219 U = cell(nGrids, 1);
220 for a = 1:nGrids
221 U{a} = cell(dim, 1);
222 I = speye(N_u{a}, N_u{a});
223 for i = 1:dim
224 E = sparse(dim,1);
225 E(i) = 1;
226 U{a}{i} = kron(I, E);
227 end
228 end
229 obj.U = U;
230
231 % Order grids
232 % u1, u2
233 Iu1 = speye(dim*N_u{1}, dim*N_u{1});
234 Iu2 = speye(dim*N_u{2}, dim*N_u{2});
235
236 Gu1 = cell2mat( {Iu1; sparse(dim*N_u{2}, dim*N_u{1})} );
237 Gu2 = cell2mat( {sparse(dim*N_u{1}, dim*N_u{2}); Iu2} );
238
239 G = {Gu1; Gu2};
240
241 %---- Grid layout -------
242 % gu1 = xp o yp;
243 % gu2 = xd o yd;
244 % gs1 = xd o yp;
245 % gs2 = xp o yd;
246 %------------------------
247
248 % Boundary restriction ops
249 e_w_u = cell(nGrids, 1);
250 e_s_u = cell(nGrids, 1);
251 e_e_u = cell(nGrids, 1);
252 e_n_u = cell(nGrids, 1);
253
254 e_w_s = cell(nGrids, 1);
255 e_s_s = cell(nGrids, 1);
256 e_e_s = cell(nGrids, 1);
257 e_n_s = cell(nGrids, 1);
258
259 e_w_u{1} = kron(ep_l{1}, Ip{2});
260 e_e_u{1} = kron(ep_r{1}, Ip{2});
261 e_s_u{1} = kron(Ip{1}, ep_l{2});
262 e_n_u{1} = kron(Ip{1}, ep_r{2});
263
264 e_w_u{2} = kron(ed_l{1}, Id{2});
265 e_e_u{2} = kron(ed_r{1}, Id{2});
266 e_s_u{2} = kron(Id{1}, ed_l{2});
267 e_n_u{2} = kron(Id{1}, ed_r{2});
268
269 e_w_s{1} = kron(ed_l{1}, Ip{2});
270 e_e_s{1} = kron(ed_r{1}, Ip{2});
271 e_s_s{1} = kron(Id{1}, ep_l{2});
272 e_n_s{1} = kron(Id{1}, ep_r{2});
273
274 e_w_s{2} = kron(ep_l{1}, Id{2});
275 e_e_s{2} = kron(ep_r{1}, Id{2});
276 e_s_s{2} = kron(Ip{1}, ed_l{2});
277 e_n_s{2} = kron(Ip{1}, ed_r{2});
278
279 % --- Metric coefficients on stress grids -------
280 x = cell(nGrids, 1);
281 y = cell(nGrids, 1);
282 J = cell(nGrids, 1);
283 x_xi = cell(nGrids, 1);
284 x_eta = cell(nGrids, 1);
285 y_xi = cell(nGrids, 1);
286 y_eta = cell(nGrids, 1);
287
288 for a = 1:nGrids
289 coords = g_u{a}.points();
290 x{a} = coords(:,1);
291 y{a} = coords(:,2);
292 end
293
294 for a = 1:nGrids
295 x_xi{a} = zeros(N_s{a}, 1);
296 y_xi{a} = zeros(N_s{a}, 1);
297 x_eta{a} = zeros(N_s{a}, 1);
298 y_eta{a} = zeros(N_s{a}, 1);
299
300 for b = 1:nGrids
301 x_xi{a} = x_xi{a} + D1_u2s{a,b}{1}*x{b};
302 y_xi{a} = y_xi{a} + D1_u2s{a,b}{1}*y{b};
303 x_eta{a} = x_eta{a} + D1_u2s{a,b}{2}*x{b};
304 y_eta{a} = y_eta{a} + D1_u2s{a,b}{2}*y{b};
305 end
306 end
307
308 for a = 1:nGrids
309 J{a} = x_xi{a}.*y_eta{a} - x_eta{a}.*y_xi{a};
310 end
311
312 K = cell(nGrids, 1);
313 for a = 1:nGrids
314 K{a} = cell(dim, dim);
315 K{a}{1,1} = y_eta{a}./J{a};
316 K{a}{1,2} = -y_xi{a}./J{a};
317 K{a}{2,1} = -x_eta{a}./J{a};
318 K{a}{2,2} = x_xi{a}./J{a};
319 end
320 % ----------------------------------------------
321
322 % --- Metric coefficients on displacement grids -------
323 x_s = cell(nGrids, 1);
324 y_s = cell(nGrids, 1);
325 J_u = cell(nGrids, 1);
326 x_xi_u = cell(nGrids, 1);
327 x_eta_u = cell(nGrids, 1);
328 y_xi_u = cell(nGrids, 1);
329 y_eta_u = cell(nGrids, 1);
330
331 for a = 1:nGrids
332 coords = g_s{a}.points();
333 x_s{a} = coords(:,1);
334 y_s{a} = coords(:,2);
335 end
336
337 for a = 1:nGrids
338 x_xi_u{a} = zeros(N_u{a}, 1);
339 y_xi_u{a} = zeros(N_u{a}, 1);
340 x_eta_u{a} = zeros(N_u{a}, 1);
341 y_eta_u{a} = zeros(N_u{a}, 1);
342
343 for b = 1:nGrids
344 x_xi_u{a} = x_xi_u{a} + D1_s2u{a,b}{1}*x_s{b};
345 y_xi_u{a} = y_xi_u{a} + D1_s2u{a,b}{1}*y_s{b};
346 x_eta_u{a} = x_eta_u{a} + D1_s2u{a,b}{2}*x_s{b};
347 y_eta_u{a} = y_eta_u{a} + D1_s2u{a,b}{2}*y_s{b};
348 end
349 end
350
351 for a = 1:nGrids
352 J_u{a} = x_xi_u{a}.*y_eta_u{a} - x_eta_u{a}.*y_xi_u{a};
353 end
354 % ----------------------------------------------
355
356 % x_u = Du*x;
357 % x_v = Dv*x;
358 % y_u = Du*y;
359 % y_v = Dv*y;
360
361 % J = x_u.*y_v - x_v.*y_u;
362
363 % K = cell(dim, dim);
364 % K{1,1} = y_v./J;
365 % K{1,2} = -y_u./J;
366 % K{2,1} = -x_v./J;
367 % K{2,2} = x_u./J;
368
369 % Physical derivatives
370 % obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
371 % obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
372
373 % Wrap around Aniosotropic Cartesian. Transformed density and stiffness
374 rho_tilde = cell(nGrids, 1);
375 PHI = cell(nGrids, 1);
376
377 for a = 1:nGrids
378 rho_tilde{a} = J_u{a}.*rho{a};
379 end
380
381 for a = 1:nGrids
382 PHI{a} = cell(dim,dim,dim,dim);
383 for i = 1:dim
384 for j = 1:dim
385 for k = 1:dim
386 for l = 1:dim
387 PHI{a}{i,j,k,l} = 0*C{a}{i,j,k,l};
388 for m = 1:dim
389 for n = 1:dim
390 PHI{a}{i,j,k,l} = PHI{a}{i,j,k,l} + J{a}.*K{a}{m,i}.*C{a}{m,j,n,l}.*K{a}{n,k};
391 end
392 end
393 end
394 end
395 end
396 end
397 end
398
399 refObj = scheme.Elastic2dStaggeredAnisotropic(g.logic, order, rho_tilde, PHI);
400
401 G = refObj.G;
402 U = refObj.U;
403 H_u = refObj.H_u;
404
405 % Volume quadrature
406 [m, n] = size(refObj.H);
407 obj.H = sparse(m, n);
408 obj.J = sparse(m, n);
409 for a = 1:nGrids
410 for i = 1:dim
411 obj.H = obj.H + G{a}*U{a}{i}*spdiag(J_u{a})*refObj.H_u{a}*U{a}{i}'*G{a}';
412 obj.J = obj.J + G{a}*U{a}{i}*spdiag(J_u{a})*U{a}{i}'*G{a}';
413 end
414 end
415 obj.Ji = inv(obj.J);
416
417 % Boundary quadratures on stress grids
418 s_w = cell(nGrids, 1);
419 s_e = cell(nGrids, 1);
420 s_s = cell(nGrids, 1);
421 s_n = cell(nGrids, 1);
422
423 % e_w_u = refObj.e_w_u;
424 % e_e_u = refObj.e_e_u;
425 % e_s_u = refObj.e_s_u;
426 % e_n_u = refObj.e_n_u;
427
428 e_w_s = refObj.e_w_s;
429 e_e_s = refObj.e_e_s;
430 e_s_s = refObj.e_s_s;
431 e_n_s = refObj.e_n_s;
432
433 for a = 1:nGrids
434 s_w{a} = sqrt((e_w_s{a}'*x_eta{a}).^2 + (e_w_s{a}'*y_eta{a}).^2);
435 s_e{a} = sqrt((e_e_s{a}'*x_eta{a}).^2 + (e_e_s{a}'*y_eta{a}).^2);
436 s_s{a} = sqrt((e_s_s{a}'*x_xi{a}).^2 + (e_s_s{a}'*y_xi{a}).^2);
437 s_n{a} = sqrt((e_n_s{a}'*x_xi{a}).^2 + (e_n_s{a}'*y_xi{a}).^2);
438 end
439
440 obj.s_w = s_w;
441 obj.s_e = s_e;
442 obj.s_s = s_s;
443 obj.s_n = s_n;
444
445 % for a = 1:nGrids
446 % obj.H_w_s{a} = refObj.H_w_s{a}*spdiag(s_w{a});
447 % obj.H_e_s{a} = refObj.H_e_s{a}*spdiag(s_e{a});
448 % obj.H_s_s{a} = refObj.H_s_s{a}*spdiag(s_s{a});
449 % obj.H_n_s{a} = refObj.H_n_s{a}*spdiag(s_n{a});
450 % end
451
452 % Restriction operators
453 obj.e_w_u = refObj.e_w_u;
454 obj.e_e_u = refObj.e_e_u;
455 obj.e_s_u = refObj.e_s_u;
456 obj.e_n_u = refObj.e_n_u;
457
458 obj.e_w_s = refObj.e_w_s;
459 obj.e_e_s = refObj.e_e_s;
460 obj.e_s_s = refObj.e_s_s;
461 obj.e_n_s = refObj.e_n_s;
462
463 % Adapt things from reference object
464 obj.D = refObj.D;
465 obj.U = refObj.U;
466 obj.G = refObj.G;
467
468 % obj.e1_w = refObj.e1_w;
469 % obj.e1_e = refObj.e1_e;
470 % obj.e1_s = refObj.e1_s;
471 % obj.e1_n = refObj.e1_n;
472
473 % obj.e2_w = refObj.e2_w;
474 % obj.e2_e = refObj.e2_e;
475 % obj.e2_s = refObj.e2_s;
476 % obj.e2_n = refObj.e2_n;
477
478 % obj.e_scalar_w = refObj.e_scalar_w;
479 % obj.e_scalar_e = refObj.e_scalar_e;
480 % obj.e_scalar_s = refObj.e_scalar_s;
481 % obj.e_scalar_n = refObj.e_scalar_n;
482
483 % e1_w = obj.e1_w;
484 % e1_e = obj.e1_e;
485 % e1_s = obj.e1_s;
486 % e1_n = obj.e1_n;
487
488 % e2_w = obj.e2_w;
489 % e2_e = obj.e2_e;
490 % e2_s = obj.e2_s;
491 % e2_n = obj.e2_n;
492
493 % obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
494 % obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
495 % obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
496 % obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
497
498 % obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
499 % obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
500 % obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
501 % obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
502
503 % obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
504 % obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
505 % obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
506 % obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
507
508 % % Physical normals
509 % e_w = obj.e_scalar_w;
510 % e_e = obj.e_scalar_e;
511 % e_s = obj.e_scalar_s;
512 % e_n = obj.e_scalar_n;
513
514 % e_w_vec = obj.e_w;
515 % e_e_vec = obj.e_e;
516 % e_s_vec = obj.e_s;
517 % e_n_vec = obj.e_n;
518
519 % nu_w = [-1,0];
520 % nu_e = [1,0];
521 % nu_s = [0,-1];
522 % nu_n = [0,1];
523
524 % obj.n_w = cell(2,1);
525 % obj.n_e = cell(2,1);
526 % obj.n_s = cell(2,1);
527 % obj.n_n = cell(2,1);
528
529 % n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
530 % n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
531 % obj.n_w{1} = spdiag(n_w_1);
532 % obj.n_w{2} = spdiag(n_w_2);
533
534 % n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
535 % n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
536 % obj.n_e{1} = spdiag(n_e_1);
537 % obj.n_e{2} = spdiag(n_e_2);
538
539 % n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
540 % n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
541 % obj.n_s{1} = spdiag(n_s_1);
542 % obj.n_s{2} = spdiag(n_s_2);
543
544 % n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
545 % n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
546 % obj.n_n{1} = spdiag(n_n_1);
547 % obj.n_n{2} = spdiag(n_n_2);
548
549 % % Operators that extract the normal component
550 % obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
551 % obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
552 % obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
553 % obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
554
555 % Misc.
556 obj.refObj = refObj;
557 obj.m = refObj.m;
558 obj.h = refObj.h;
559 obj.order = order;
560 obj.grid = g;
561 obj.dim = dim;
562 obj.nGrids = nGrids;
563
564 end
565
566
567 % Closure functions return the operators applied to the own domain to close the boundary
568 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
569 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
570 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
571 % on the first component. Can also be e.g.
572 % {'normal', 'd'} or {'tangential', 't'} for conditions on
573 % tangential/normal component.
574 % data is a function returning the data that should be applied at the boundary.
575 % neighbour_scheme is an instance of Scheme that should be interfaced to.
576 % neighbour_boundary is a string specifying which boundary to interface to.
577
578 % For displacement bc:
579 % bc = {comp, 'd', dComps},
580 % where
581 % dComps = vector of components with displacement BC. Default: 1:dim.
582 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
583 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
584 default_arg('tuning', 1.0);
585 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
586
587 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
588
589 type = bc{2};
590
591 switch type
592 case {'F','f','Free','free','traction','Traction','t','T'}
593 s = obj.(['s_' boundary]);
594 s = spdiag(cell2mat(s));
595 penalty = penalty*s;
596 end
597 end
598
599 % type Struct that specifies the interface coupling.
600 % Fields:
601 % -- tuning: penalty strength, defaults to 1.0
602 % -- interpolation: type of interpolation, default 'none'
603 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
604
605 defaultType.tuning = 1.0;
606 defaultType.interpolation = 'none';
607 default_struct('type', defaultType);
608
609 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
610 end
611
612 % Returns h11 for the boundary specified by the string boundary.
613 % op -- string
614 function h11 = getBorrowing(obj, boundary)
615 assertIsMember(boundary, {'w', 'e', 's', 'n'})
616
617 switch boundary
618 case {'w','e'}
619 h11 = obj.refObj.h11{1};
620 case {'s', 'n'}
621 h11 = obj.refObj.h11{2};
622 end
623 end
624
625 % Returns the outward unit normal vector for the boundary specified by the string boundary.
626 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
627 function n = getNormal(obj, boundary)
628 error('Not implemented');
629 assertIsMember(boundary, {'w', 'e', 's', 'n'})
630
631 n = obj.(['n_' boundary]);
632 end
633
634 % Returns the boundary operator op for the boundary specified by the string boundary.
635 % op -- string
636 function o = getBoundaryOperator(obj, op, boundary)
637 assertIsMember(boundary, {'w', 'e', 's', 'n'})
638 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
639
640 o = obj.([op, '_', boundary]);
641
642 end
643
644 % Returns the boundary operator op for the boundary specified by the string boundary.
645 % op -- string
646 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
647 assertIsMember(boundary, {'w', 'e', 's', 'n'})
648 assertIsMember(op, {'e_u', 'e_s'})
649
650 switch op
651 case 'e_u'
652 o = obj.(['e_', boundary, '_u']);
653 case 'e_s'
654 o = obj.(['e_', boundary, '_s']);
655 end
656
657 end
658
659 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
660 % Formula: tau_i = T_ij u_j
661 % op -- string
662 function T = getBoundaryTractionOperator(obj, boundary)
663 assertIsMember(boundary, {'w', 'e', 's', 'n'})
664
665 T = obj.(['T', '_', boundary]);
666 end
667
668 % Returns square boundary quadrature matrix, of dimension
669 % corresponding to the number of boundary unknowns
670 %
671 % boundary -- string
672 function H = getBoundaryQuadrature(obj, boundary)
673 assertIsMember(boundary, {'w', 'e', 's', 'n'})
674
675 H = obj.getBoundaryQuadratureForScalarField(boundary);
676 I_dim = speye(obj.dim, obj.dim);
677 H = kron(H, I_dim);
678 end
679
680 % Returns square boundary quadrature matrix, of dimension
681 % corresponding to the number of boundary grid points
682 %
683 % boundary -- string
684 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
685 assertIsMember(boundary, {'w', 'e', 's', 'n'})
686
687 H_b = obj.(['H_', boundary]);
688 end
689
690 function N = size(obj)
691 N = length(obj.D);
692 end
693 end
694 end