Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1331:60c875c18de3 feature/D2_boundary_opt

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with feature/poroelastic for Elastic schemes
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Thu, 10 Mar 2022 16:54:26 +0100
parents 412b8ceafbc6
children df8c71b80c33
comparison
equal deleted inserted replaced
1330:855871e0b852 1331:60c875c18de3
1 classdef Elastic2dCurvilinearAnisotropic < scheme.Scheme
2
3 % Discretizes the elastic wave equation:
4 % rho u_{i,tt} = dj C_{ijkl} dk u_j
5 % in curvilinear coordinates.
6 % opSet should be cell array of opSets, one per dimension. This
7 % is useful if we have periodic BC in one direction.
8 % Assumes fully compatible operators.
9
10 properties
11 m % Number of points in each direction, possibly a vector
12 h % Grid spacing
13
14 grid
15 dim
16
17 order % Order of accuracy for the approximation
18
19 % Diagonal matrices for variable coefficients
20 J, Ji
21 RHO % Density
22 C % Elastic stiffness tensor
23
24 D % Total operator
25
26 K % Transformation gradient
27 Dx, Dy % Physical derivatives
28 sigma % Cell matrix of physical stress operators
29 n_w, n_e, n_s, n_n % Physical normals
30 tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
31
32 % Boundary operators in cell format, used for BC
33 T_w, T_e, T_s, T_n
34
35 % Traction operators
36 tau_w, tau_e, tau_s, tau_n % Return vector field
37 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
38 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
39 tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
40 tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
41
42 % Inner products
43 H
44
45 % Boundary inner products (for scalar field)
46 H_w, H_e, H_s, H_n
47
48 % Surface Jacobian vectors
49 s_w, s_e, s_s, s_n
50
51 % Boundary restriction operators
52 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
53 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
54 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
55 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
56 en_w, en_e, en_s, en_n % Act on vector field, return normal component
57 et_w, et_e, et_s, et_n % Act on vector field, return tangential component
58
59 % E{i}^T picks out component i
60 E
61
62 % Elastic2dVariableAnisotropic object for reference domain
63 refObj
64 end
65
66 methods
67
68 % The coefficients can either be function handles or grid functions
69 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
70 function obj = Elastic2dCurvilinearAnisotropic(g, order, rho, C, opSet, optFlag, hollow)
71 default_arg('hollow', false);
72 default_arg('rho', @(x,y) 0*x+1);
73 default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
74 default_arg('optFlag', false);
75 dim = 2;
76
77 C_default = cell(dim,dim,dim,dim);
78 for i = 1:dim
79 for j = 1:dim
80 for k = 1:dim
81 for l = 1:dim
82 C_default{i,j,k,l} = @(x,y) 0*x ;
83 end
84 end
85 end
86 end
87 default_arg('C', C_default);
88
89 assert(isa(g, 'grid.Curvilinear'));
90
91 if isa(rho, 'function_handle')
92 rho = grid.evalOn(g, rho);
93 end
94
95 C_mat = cell(dim,dim,dim,dim);
96 for i = 1:dim
97 for j = 1:dim
98 for k = 1:dim
99 for l = 1:dim
100 if isa(C{i,j,k,l}, 'function_handle')
101 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
102 end
103 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
104 end
105 end
106 end
107 end
108 obj.C = C_mat;
109
110 m = g.size();
111 m_tot = g.N();
112
113 % 1D operators
114 m_u = m(1);
115 m_v = m(2);
116 ops_u = opSet{1}(m_u, {0, 1}, order);
117 ops_v = opSet{2}(m_v, {0, 1}, order);
118
119 h_u = ops_u.h;
120 h_v = ops_v.h;
121
122 I_u = speye(m_u);
123 I_v = speye(m_v);
124
125 D1_u = ops_u.D1;
126 H_u = ops_u.H;
127 Hi_u = ops_u.HI;
128 e_l_u = ops_u.e_l;
129 e_r_u = ops_u.e_r;
130 d1_l_u = ops_u.d1_l;
131 d1_r_u = ops_u.d1_r;
132
133 D1_v = ops_v.D1;
134 H_v = ops_v.H;
135 Hi_v = ops_v.HI;
136 e_l_v = ops_v.e_l;
137 e_r_v = ops_v.e_r;
138 d1_l_v = ops_v.d1_l;
139 d1_r_v = ops_v.d1_r;
140
141
142 % Logical operators
143 Du = kr(D1_u,I_v);
144 Dv = kr(I_u,D1_v);
145
146 e_w = kr(e_l_u,I_v);
147 e_e = kr(e_r_u,I_v);
148 e_s = kr(I_u,e_l_v);
149 e_n = kr(I_u,e_r_v);
150
151 % Metric coefficients
152 coords = g.points();
153 x = coords(:,1);
154 y = coords(:,2);
155
156 x_u = Du*x;
157 x_v = Dv*x;
158 y_u = Du*y;
159 y_v = Dv*y;
160
161 J = x_u.*y_v - x_v.*y_u;
162
163 K = cell(dim, dim);
164 K{1,1} = y_v./J;
165 K{1,2} = -y_u./J;
166 K{2,1} = -x_v./J;
167 K{2,2} = x_u./J;
168 obj.K = K;
169
170 % Physical derivatives
171 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
172 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
173
174 % Wrap around Aniosotropic Cartesian
175 rho_tilde = J.*rho;
176
177 PHI = cell(dim,dim,dim,dim);
178 for i = 1:dim
179 for j = 1:dim
180 for k = 1:dim
181 for l = 1:dim
182 PHI{i,j,k,l} = 0*C{i,j,k,l};
183 for m = 1:dim
184 for n = 1:dim
185 PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,i}.*C{m,j,n,l}.*K{n,k};
186 end
187 end
188 end
189 end
190 end
191 end
192
193 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
194 refObj = scheme.Elastic2dVariableAnisotropic(gRef, order, rho_tilde, PHI, opSet, [], hollow);
195
196 %---- Set object properties ------
197 obj.RHO = spdiag(rho);
198
199 % Volume quadrature
200 obj.J = spdiag(J);
201 obj.Ji = spdiag(1./J);
202 obj.H = obj.J*kr(H_u,H_v);
203
204 % Boundary quadratures
205 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
206 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
207 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
208 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
209 obj.s_w = s_w;
210 obj.s_e = s_e;
211 obj.s_s = s_s;
212 obj.s_n = s_n;
213
214 obj.H_w = H_v*spdiag(s_w);
215 obj.H_e = H_v*spdiag(s_e);
216 obj.H_s = H_u*spdiag(s_s);
217 obj.H_n = H_u*spdiag(s_n);
218
219 % Restriction operators
220 obj.e_w = refObj.e_w;
221 obj.e_e = refObj.e_e;
222 obj.e_s = refObj.e_s;
223 obj.e_n = refObj.e_n;
224
225 % Adapt things from reference object
226 obj.D = refObj.D;
227 obj.E = refObj.E;
228
229 obj.e1_w = refObj.e1_w;
230 obj.e1_e = refObj.e1_e;
231 obj.e1_s = refObj.e1_s;
232 obj.e1_n = refObj.e1_n;
233
234 obj.e2_w = refObj.e2_w;
235 obj.e2_e = refObj.e2_e;
236 obj.e2_s = refObj.e2_s;
237 obj.e2_n = refObj.e2_n;
238
239 obj.e_scalar_w = refObj.e_scalar_w;
240 obj.e_scalar_e = refObj.e_scalar_e;
241 obj.e_scalar_s = refObj.e_scalar_s;
242 obj.e_scalar_n = refObj.e_scalar_n;
243
244 e1_w = obj.e1_w;
245 e1_e = obj.e1_e;
246 e1_s = obj.e1_s;
247 e1_n = obj.e1_n;
248
249 e2_w = obj.e2_w;
250 e2_e = obj.e2_e;
251 e2_s = obj.e2_s;
252 e2_n = obj.e2_n;
253
254 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
255 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
256 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
257 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
258
259 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
260 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
261 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
262 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
263
264 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
265 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
266 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
267 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
268
269 % Physical normals
270 e_w = obj.e_scalar_w;
271 e_e = obj.e_scalar_e;
272 e_s = obj.e_scalar_s;
273 e_n = obj.e_scalar_n;
274
275 e_w_vec = obj.e_w;
276 e_e_vec = obj.e_e;
277 e_s_vec = obj.e_s;
278 e_n_vec = obj.e_n;
279
280 nu_w = [-1,0];
281 nu_e = [1,0];
282 nu_s = [0,-1];
283 nu_n = [0,1];
284
285 obj.n_w = cell(2,1);
286 obj.n_e = cell(2,1);
287 obj.n_s = cell(2,1);
288 obj.n_n = cell(2,1);
289
290 % Compute normal and rotate (exactly!) 90 degrees counter-clockwise to get tangent
291 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
292 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
293 obj.n_w{1} = spdiag(n_w_1);
294 obj.n_w{2} = spdiag(n_w_2);
295 obj.tangent_w = {-obj.n_w{2}, obj.n_w{1}};
296
297 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
298 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
299 obj.n_e{1} = spdiag(n_e_1);
300 obj.n_e{2} = spdiag(n_e_2);
301 obj.tangent_e = {-obj.n_e{2}, obj.n_e{1}};
302
303 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
304 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
305 obj.n_s{1} = spdiag(n_s_1);
306 obj.n_s{2} = spdiag(n_s_2);
307 obj.tangent_s = {-obj.n_s{2}, obj.n_s{1}};
308
309 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
310 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
311 obj.n_n{1} = spdiag(n_n_1);
312 obj.n_n{2} = spdiag(n_n_2);
313 obj.tangent_n = {-obj.n_n{2}, obj.n_n{1}};
314
315 % Operators that extract the normal component
316 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
317 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
318 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
319 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
320
321 % Operators that extract the tangential component
322 obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')';
323 obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')';
324 obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')';
325 obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')';
326
327 obj.tau_n_w = (obj.n_w{1}*obj.tau1_w' + obj.n_w{2}*obj.tau2_w')';
328 obj.tau_n_e = (obj.n_e{1}*obj.tau1_e' + obj.n_e{2}*obj.tau2_e')';
329 obj.tau_n_s = (obj.n_s{1}*obj.tau1_s' + obj.n_s{2}*obj.tau2_s')';
330 obj.tau_n_n = (obj.n_n{1}*obj.tau1_n' + obj.n_n{2}*obj.tau2_n')';
331
332 obj.tau_t_w = (obj.tangent_w{1}*obj.tau1_w' + obj.tangent_w{2}*obj.tau2_w')';
333 obj.tau_t_e = (obj.tangent_e{1}*obj.tau1_e' + obj.tangent_e{2}*obj.tau2_e')';
334 obj.tau_t_s = (obj.tangent_s{1}*obj.tau1_s' + obj.tangent_s{2}*obj.tau2_s')';
335 obj.tau_t_n = (obj.tangent_n{1}*obj.tau1_n' + obj.tangent_n{2}*obj.tau2_n')';
336
337 % Stress operators
338 sigma = cell(dim, dim);
339 D1 = {obj.Dx, obj.Dy};
340 E = obj.E;
341 N = length(obj.RHO);
342 for i = 1:dim
343 for j = 1:dim
344 sigma{i,j} = sparse(N,2*N);
345 for k = 1:dim
346 for l = 1:dim
347 sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}';
348 end
349 end
350 end
351 end
352 obj.sigma = sigma;
353
354 % Misc.
355 obj.refObj = refObj;
356 obj.m = refObj.m;
357 obj.h = refObj.h;
358 obj.order = order;
359 obj.grid = g;
360 obj.dim = dim;
361
362 end
363
364
365 % Closure functions return the operators applied to the own domain to close the boundary
366 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
367 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
368 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
369 % on the first component. Can also be e.g.
370 % {'normal', 'd'} or {'tangential', 't'} for conditions on
371 % tangential/normal component.
372 % data is a function returning the data that should be applied at the boundary.
373 % neighbour_scheme is an instance of Scheme that should be interfaced to.
374 % neighbour_boundary is a string specifying which boundary to interface to.
375
376 % For displacement bc:
377 % bc = {comp, 'd', dComps},
378 % where
379 % dComps = vector of components with displacement BC. Default: 1:dim.
380 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
381 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
382 default_arg('tuning', 1.0);
383 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
384
385 component = bc{1};
386 type = bc{2};
387
388 switch component
389
390 % If conditions on Cartesian components
391 case {1,2}
392 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
393
394 % If conditions on normal or tangential components
395 case {'n', 't'}
396
397 switch component
398 case 'n'
399 en = obj.getBoundaryOperator('en', boundary);
400 case 't'
401 en = obj.getBoundaryOperator('et', boundary);
402 end
403 e1 = obj.getBoundaryOperator('e1', boundary);
404
405 bc1 = {1, type};
406 [c1, p1] = obj.refObj.boundary_condition(boundary, bc1, tuning);
407 bc2 = {2, type};
408 c2 = obj.refObj.boundary_condition(boundary, bc2, tuning);
409
410 switch type
411 case {'F','f','Free','free','traction','Traction','t','T'}
412 closure = en*en'*(c1+c2);
413 penalty = en*e1'*p1;
414 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
415 [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
416 end
417
418 end
419
420 switch type
421 case {'F','f','Free','free','traction','Traction','t','T'}
422
423 s = obj.(['s_' boundary]);
424 s = spdiag(s);
425 penalty = penalty*s;
426
427 end
428 end
429
430 function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
431 u = obj;
432
433 component = bc{1};
434 type = bc{2};
435
436 switch component
437 case 'n'
438 en = u.getBoundaryOperator('en', boundary);
439 tau_n = u.getBoundaryOperator('tau_n', boundary);
440 N = u.getNormal(boundary);
441 case 't'
442 en = u.getBoundaryOperator('et', boundary);
443 tau_n = u.getBoundaryOperator('tau_t', boundary);
444 N = u.getTangent(boundary);
445 end
446
447 % Operators
448 e = u.getBoundaryOperatorForScalarField('e', boundary);
449 h11 = u.getBorrowing(boundary);
450 n = u.getNormal(boundary);
451
452 C = u.C;
453 Ji = u.Ji;
454 s = spdiag(u.(['s_' boundary]));
455 m_tot = u.grid.N();
456
457 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
458 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
459
460 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
461 dim = u.dim;
462
463 % Preallocate
464 [~, m_int] = size(H_gamma);
465 closure = sparse(dim*m_tot, dim*m_tot);
466 penalty = sparse(dim*m_tot, m_int);
467
468 % Term 1: The symmetric term
469 Z = sparse(m_int, m_int);
470 for i = 1:dim
471 for j = 1:dim
472 for l = 1:dim
473 for k = 1:dim
474 Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
475 end
476 end
477 end
478 end
479
480 Z = -tuning*dim*1/h11*s*Z;
481 closure = closure + en*H_gamma*Z*en';
482 penalty = penalty - en*H_gamma*Z;
483
484 % Term 2: The symmetrizing term
485 closure = closure + tau_n*H_gamma*en';
486 penalty = penalty - tau_n*H_gamma;
487
488 % Multiply all terms by inverse of density x quadraure
489 closure = RHOi*Hi*closure;
490 penalty = RHOi*Hi*penalty;
491 end
492
493 % type Struct that specifies the interface coupling.
494 % Fields:
495 % -- tuning: penalty strength, defaults to 1.0
496 % -- interpolation: type of interpolation, default 'none'
497 function [closure, penalty, forcingPenalties] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
498
499 defaultType.tuning = 1.0;
500 defaultType.interpolation = 'none';
501 defaultType.type = 'standard';
502 default_struct('type', defaultType);
503
504 forcingPenalties = [];
505
506 switch type.type
507 case 'standard'
508 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
509 case 'normalTangential'
510 [closure, penalty, forcingPenalties] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
511 case 'frictionalFault'
512 [closure, penalty, forcingPenalties] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
513 end
514
515 end
516
517 function [closure, penalty, forcingPenalties] = interfaceFrictionalFault(obj,boundary,neighbour_scheme,neighbour_boundary,type)
518 tuning = type.tuning;
519
520 % u denotes the solution in the own domain
521 % v denotes the solution in the neighbour domain
522
523 forcingPenalties = cell(1, 1);
524 u = obj;
525 v = neighbour_scheme;
526
527 % Operators, u side
528 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
529 en_u = u.getBoundaryOperator('en', boundary);
530 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
531 h11_u = u.getBorrowing(boundary);
532 n_u = u.getNormal(boundary);
533
534 C_u = u.C;
535 Ji_u = u.Ji;
536 s_u = spdiag(u.(['s_' boundary]));
537 m_tot_u = u.grid.N();
538
539 % Operators, v side
540 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
541 en_v = v.getBoundaryOperator('en', neighbour_boundary);
542 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
543 h11_v = v.getBorrowing(neighbour_boundary);
544 n_v = v.getNormal(neighbour_boundary);
545
546 C_v = v.C;
547 Ji_v = v.Ji;
548 s_v = spdiag(v.(['s_' neighbour_boundary]));
549 m_tot_v = v.grid.N();
550
551 % Operators that are only required for own domain
552 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
553 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
554
555 % Shared operators
556 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
557 dim = u.dim;
558
559 % Preallocate
560 [~, m_int] = size(H_gamma);
561 closure = sparse(dim*m_tot_u, dim*m_tot_u);
562 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
563
564 % Continuity of normal displacement, term 1: The symmetric term
565 Z_u = sparse(m_int, m_int);
566 Z_v = sparse(m_int, m_int);
567 for i = 1:dim
568 for j = 1:dim
569 for l = 1:dim
570 for k = 1:dim
571 Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
572 Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
573 end
574 end
575 end
576 end
577
578 Z = -tuning*dim*( 1/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
579 closure = closure + en_u*H_gamma*Z*en_u';
580 penalty = penalty + en_u*H_gamma*Z*en_v';
581
582 % Continuity of normal displacement, term 2: The symmetrizing term
583 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
584 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
585
586 % Continuity of normal traction
587 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
588 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
589 forcing_tau_n = 1/2*en_u*H_gamma;
590
591 % Multiply all normal component terms by inverse of density x quadraure
592 closure = RHOi*Hi*closure;
593 penalty = RHOi*Hi*penalty;
594 forcing_tau_n = RHOi*Hi*forcing_tau_n;
595
596 % ---- Tangential tractions are imposed just like traction BC ------
597 closure = closure + obj.boundary_condition(boundary, {'t', 't'});
598
599 forcingPenalties{1} = forcing_tau_n;
600
601 end
602
603 % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential
604 % coordinate system. For the isotropic case, the components decouple nicely.
605 % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
606 function [closure, penalty, forcingPenalties, Zt] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
607 tuning = type.tuning;
608
609 % u denotes the solution in the own domain
610 % v denotes the solution in the neighbour domain
611
612 forcingPenalties = cell(2, 1);
613 u = obj;
614 v = neighbour_scheme;
615
616 % Operators, u side
617 e_u = u.getBoundaryOperatorForScalarField('e', boundary);
618 en_u = u.getBoundaryOperator('en', boundary);
619 et_u = u.getBoundaryOperator('et', boundary);
620 tau_n_u = u.getBoundaryOperator('tau_n', boundary);
621 tau_t_u = u.getBoundaryOperator('tau_t', boundary);
622 h11_u = u.getBorrowing(boundary);
623 n_u = u.getNormal(boundary);
624 t_u = u.getTangent(boundary);
625
626 C_u = u.C;
627 Ji_u = u.Ji;
628 s_u = spdiag(u.(['s_' boundary]));
629 m_tot_u = u.grid.N();
630
631 % Operators, v side
632 e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
633 en_v = v.getBoundaryOperator('en', neighbour_boundary);
634 et_v = v.getBoundaryOperator('et', neighbour_boundary);
635 tau_n_v = v.getBoundaryOperator('tau_n', neighbour_boundary);
636 tau_t_v = v.getBoundaryOperator('tau_t', neighbour_boundary);
637 h11_v = v.getBorrowing(neighbour_boundary);
638 n_v = v.getNormal(neighbour_boundary);
639 t_v = v.getTangent(neighbour_boundary);
640
641 C_v = v.C;
642 Ji_v = v.Ji;
643 s_v = spdiag(v.(['s_' neighbour_boundary]));
644 m_tot_v = v.grid.N();
645
646 % Operators that are only required for own domain
647 Hi = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
648 RHOi = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
649
650 % Shared operators
651 H_gamma = u.getBoundaryQuadratureForScalarField(boundary);
652 dim = u.dim;
653
654 % Preallocate
655 [~, m_int] = size(H_gamma);
656 closure = sparse(dim*m_tot_u, dim*m_tot_u);
657 penalty = sparse(dim*m_tot_u, dim*m_tot_v);
658
659 % -- Continuity of displacement, term 1: The symmetric term ---
660 Zn_u = sparse(m_int, m_int);
661 Zn_v = sparse(m_int, m_int);
662 Zt_u = sparse(m_int, m_int);
663 Zt_v = sparse(m_int, m_int);
664 for i = 1:dim
665 for j = 1:dim
666 for l = 1:dim
667 for k = 1:dim
668 % Penalty strength for normal component
669 Zn_u = Zn_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
670 Zn_v = Zn_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
671
672 % Penalty strength for tangential component
673 Zt_u = Zt_u + n_u{i}*t_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*t_u{l};
674 Zt_v = Zt_v + n_v{i}*t_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*t_v{l};
675 end
676 end
677 end
678 end
679
680 Zn = -tuning*dim*( 1/(4*h11_u)*s_u*Zn_u + 1/(4*h11_v)*s_v*Zn_v );
681 Zt = -tuning*dim*( 1/(4*h11_u)*s_u*Zt_u + 1/(4*h11_v)*s_v*Zt_v );
682
683 % Continuity of normal component
684 closure = closure + en_u*H_gamma*Zn*en_u';
685 penalty = penalty + en_u*H_gamma*Zn*en_v';
686 forcing_u_n = -en_u*H_gamma*Zn;
687
688 % Continuity of tangential component
689 closure = closure + et_u*H_gamma*Zt*et_u';
690 penalty = penalty + et_u*H_gamma*Zt*et_v';
691 forcing_u_t = -et_u*H_gamma*Zt;
692 %------------------------------------------------------------------
693
694 % --- Continuity of displacement, term 2: The symmetrizing term
695
696 % Continuity of normal displacement
697 closure = closure + 1/2*tau_n_u*H_gamma*en_u';
698 penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
699 forcing_u_n = forcing_u_n - 1/2*tau_n_u*H_gamma;
700
701 % Continuity of tangential displacement
702 closure = closure + 1/2*tau_t_u*H_gamma*et_u';
703 penalty = penalty + 1/2*tau_t_u*H_gamma*et_v';
704 forcing_u_t = forcing_u_t - 1/2*tau_t_u*H_gamma;
705 % ------------------------------------------------------------------
706
707 % --- Continuity of tractions -----------------------------
708
709 % Continuity of normal traction
710 closure = closure - 1/2*en_u*H_gamma*tau_n_u';
711 penalty = penalty + 1/2*en_u*H_gamma*tau_n_v';
712
713 % Continuity of tangential traction
714 closure = closure - 1/2*et_u*H_gamma*tau_t_u';
715 penalty = penalty + 1/2*et_u*H_gamma*tau_t_v';
716 %--------------------------------------------------------------------
717
718 % Multiply all terms by inverse of density x quadraure
719 closure = RHOi*Hi*closure;
720 penalty = RHOi*Hi*penalty;
721 forcing_u_n = RHOi*Hi*forcing_u_n;
722 forcing_u_t = RHOi*Hi*forcing_u_t;
723
724 forcingPenalties{1} = forcing_u_n;
725 forcingPenalties{2} = forcing_u_t;
726
727 end
728
729
730 % Returns h11 for the boundary specified by the string boundary.
731 % op -- string
732 function h11 = getBorrowing(obj, boundary)
733 assertIsMember(boundary, {'w', 'e', 's', 'n'})
734
735 switch boundary
736 case {'w','e'}
737 h11 = obj.refObj.h11{1};
738 case {'s', 'n'}
739 h11 = obj.refObj.h11{2};
740 end
741 end
742
743 % Returns the outward unit normal vector for the boundary specified by the string boundary.
744 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
745 function n = getNormal(obj, boundary)
746 assertIsMember(boundary, {'w', 'e', 's', 'n'})
747
748 n = obj.(['n_' boundary]);
749 end
750
751 % Returns the unit tangent vector for the boundary specified by the string boundary.
752 % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
753 function t = getTangent(obj, boundary)
754 assertIsMember(boundary, {'w', 'e', 's', 'n'})
755
756 t = obj.(['tangent_' boundary]);
757 end
758
759 % Returns the boundary operator op for the boundary specified by the string boundary.
760 % op -- string
761 function o = getBoundaryOperator(obj, op, boundary)
762 assertIsMember(boundary, {'w', 'e', 's', 'n'})
763 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
764
765 o = obj.([op, '_', boundary]);
766
767 end
768
769 % Returns the boundary operator op for the boundary specified by the string boundary.
770 % op -- string
771 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
772 assertIsMember(boundary, {'w', 'e', 's', 'n'})
773 assertIsMember(op, {'e'})
774
775 switch op
776
777 case 'e'
778 o = obj.(['e_scalar', '_', boundary]);
779 end
780
781 end
782
783 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
784 % Formula: tau_i = T_ij u_j
785 % op -- string
786 function T = getBoundaryTractionOperator(obj, boundary)
787 assertIsMember(boundary, {'w', 'e', 's', 'n'})
788
789 T = obj.(['T', '_', boundary]);
790 end
791
792 % Returns square boundary quadrature matrix, of dimension
793 % corresponding to the number of boundary unknowns
794 %
795 % boundary -- string
796 function H = getBoundaryQuadrature(obj, boundary)
797 assertIsMember(boundary, {'w', 'e', 's', 'n'})
798
799 H = obj.getBoundaryQuadratureForScalarField(boundary);
800 I_dim = speye(obj.dim, obj.dim);
801 H = kron(H, I_dim);
802 end
803
804 % Returns square boundary quadrature matrix, of dimension
805 % corresponding to the number of boundary grid points
806 %
807 % boundary -- string
808 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
809 assertIsMember(boundary, {'w', 'e', 's', 'n'})
810
811 H_b = obj.(['H_', boundary]);
812 end
813
814 function N = size(obj)
815 N = obj.dim*prod(obj.m);
816 end
817 end
818 end