Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.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 19d821ddc108 70939ea9a71f
children
comparison
equal deleted inserted replaced
1330:855871e0b852 1331:60c875c18de3
12 grid 12 grid
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 variable 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.
53 end 59 end
54 60
55 methods 61 methods
56 62
57 % The coefficients can either be function handles or grid functions 63 % The coefficients can either be function handles or grid functions
58 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet) 64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
65 function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet, optFlag)
59 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); 66 default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
60 default_arg('lambda', @(x,y) 0*x+1); 67 default_arg('lambda', @(x,y) 0*x+1);
61 default_arg('mu', @(x,y) 0*x+1); 68 default_arg('mu', @(x,y) 0*x+1);
62 default_arg('rho', @(x,y) 0*x+1); 69 default_arg('rho', @(x,y) 0*x+1);
70 default_arg('optFlag', false);
63 dim = 2; 71 dim = 2;
64 72
65 assert(isa(g, 'grid.Cartesian')) 73 assert(isa(g, 'grid.Cartesian'))
66 74
67 if isa(lambda, 'function_handle') 75 if isa(lambda, 'function_handle')
103 I = cell(dim,1); 111 I = cell(dim,1);
104 D1 = cell(dim,1); 112 D1 = cell(dim,1);
105 D2 = cell(dim,1); 113 D2 = cell(dim,1);
106 H = cell(dim,1); 114 H = cell(dim,1);
107 Hi = cell(dim,1); 115 Hi = cell(dim,1);
108 e_l = cell(dim,1); 116 e_0 = cell(dim,1);
109 e_r = cell(dim,1); 117 e_m = cell(dim,1);
110 d1_l = cell(dim,1); 118 d1_0 = cell(dim,1);
111 d1_r = cell(dim,1); 119 d1_m = cell(dim,1);
112 120
113 for i = 1:dim 121 for i = 1:dim
114 I{i} = speye(m(i)); 122 I{i} = speye(m(i));
115 D1{i} = ops{i}.D1; 123 D1{i} = ops{i}.D1;
116 D2{i} = ops{i}.D2; 124 D2{i} = ops{i}.D2;
117 H{i} = ops{i}.H; 125 H{i} = ops{i}.H;
118 Hi{i} = ops{i}.HI; 126 Hi{i} = ops{i}.HI;
119 e_l{i} = ops{i}.e_l; 127 e_0{i} = ops{i}.e_l;
120 e_r{i} = ops{i}.e_r; 128 e_m{i} = ops{i}.e_r;
121 d1_l{i} = ops{i}.d1_l; 129 d1_0{i} = ops{i}.d1_l;
122 d1_r{i} = ops{i}.d1_r; 130 d1_m{i} = ops{i}.d1_r;
123 end 131 end
124 132
125 %====== Assemble full operators ======== 133 %====== Assemble full operators ========
126 LAMBDA = spdiag(lambda); 134 LAMBDA = spdiag(lambda);
127 obj.LAMBDA = LAMBDA; 135 obj.LAMBDA = LAMBDA;
132 obj.RHOi = inv(RHO); 140 obj.RHOi = inv(RHO);
133 141
134 obj.D1 = cell(dim,1); 142 obj.D1 = cell(dim,1);
135 obj.D2_lambda = cell(dim,1); 143 obj.D2_lambda = cell(dim,1);
136 obj.D2_mu = cell(dim,1); 144 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 145
142 % D1 146 % D1
143 obj.D1{1} = kron(D1{1},I{2}); 147 obj.D1{1} = kron(D1{1},I{2});
144 obj.D1{2} = kron(I{1},D1{2}); 148 obj.D1{2} = kron(I{1},D1{2});
145 149
146 % Boundary operators 150 % Boundary restriction operators
147 obj.e_l{1} = kron(e_l{1},I{2}); 151 e_l = cell(dim,1);
148 obj.e_l{2} = kron(I{1},e_l{2}); 152 e_r = cell(dim,1);
149 obj.e_r{1} = kron(e_r{1},I{2}); 153 e_l{1} = kron(e_0{1}, I{2});
150 obj.e_r{2} = kron(I{1},e_r{2}); 154 e_l{2} = kron(I{1}, e_0{2});
151 155 e_r{1} = kron(e_m{1}, I{2});
152 obj.d1_l{1} = kron(d1_l{1},I{2}); 156 e_r{2} = kron(I{1}, e_m{2});
153 obj.d1_l{2} = kron(I{1},d1_l{2}); 157
154 obj.d1_r{1} = kron(d1_r{1},I{2}); 158 e_scalar_w = e_l{1};
155 obj.d1_r{2} = kron(I{1},d1_r{2}); 159 e_scalar_e = e_r{1};
156 160 e_scalar_s = e_l{2};
157 % D2 161 e_scalar_n = e_r{2};
158 for i = 1:dim 162
159 obj.D2_lambda{i} = sparse(m_tot); 163 I_dim = speye(dim, dim);
160 obj.D2_mu{i} = sparse(m_tot); 164 e_w = kron(e_scalar_w, I_dim);
161 end 165 e_e = kron(e_scalar_e, I_dim);
162 ind = grid.funcToMatrix(g, 1:m_tot); 166 e_s = kron(e_scalar_s, I_dim);
163 167 e_n = kron(e_scalar_n, I_dim);
164 for i = 1:m(2) 168
165 D_lambda = D2{1}(lambda(ind(:,i))); 169 % Boundary derivatives
166 D_mu = D2{1}(mu(ind(:,i))); 170 d1_l = cell(dim,1);
167 171 d1_r = cell(dim,1);
168 p = ind(:,i); 172 d1_l{1} = kron(d1_0{1}, I{2});
169 obj.D2_lambda{1}(p,p) = D_lambda; 173 d1_l{2} = kron(I{1}, d1_0{2});
170 obj.D2_mu{1}(p,p) = D_mu; 174 d1_r{1} = kron(d1_m{1}, I{2});
171 end 175 d1_r{2} = kron(I{1}, d1_m{2});
172 176
173 for i = 1:m(1)
174 D_lambda = D2{2}(lambda(ind(i,:)));
175 D_mu = D2{2}(mu(ind(i,:)));
176
177 p = ind(i,:);
178 obj.D2_lambda{2}(p,p) = D_lambda;
179 obj.D2_mu{2}(p,p) = D_mu;
180 end
181
182 % Quadratures
183 obj.H = kron(H{1},H{2});
184 obj.Hi = inv(obj.H);
185 obj.H_boundary = cell(dim,1);
186 obj.H_boundary{1} = H{2};
187 obj.H_boundary{2} = H{1};
188 obj.H_1D = {H{1}, H{2}};
189 177
190 % E{i}^T picks out component i. 178 % E{i}^T picks out component i.
191 E = cell(dim,1); 179 E = cell(dim,1);
192 I = speye(m_tot,m_tot); 180 I = speye(m_tot,m_tot);
193 for i = 1:dim 181 for i = 1:dim
194 e = sparse(dim,1); 182 e = sparse(dim,1);
195 e(i) = 1; 183 e(i) = 1;
196 E{i} = kron(I,e); 184 E{i} = kron(I,e);
197 end 185 end
198 obj.E = E; 186 obj.E = E;
187
188 e1_w = (e_scalar_w'*E{1}')';
189 e1_e = (e_scalar_e'*E{1}')';
190 e1_s = (e_scalar_s'*E{1}')';
191 e1_n = (e_scalar_n'*E{1}')';
192
193 e2_w = (e_scalar_w'*E{2}')';
194 e2_e = (e_scalar_e'*E{2}')';
195 e2_s = (e_scalar_s'*E{2}')';
196 e2_n = (e_scalar_n'*E{2}')';
197
198
199 % D2
200 for i = 1:dim
201 obj.D2_lambda{i} = sparse(m_tot, m_tot);
202 obj.D2_mu{i} = sparse(m_tot, m_tot);
203 end
204 ind = grid.funcToMatrix(g, 1:m_tot);
205
206 for i = 1:m(2)
207 D_lambda = D2{1}(lambda(ind(:,i)));
208 D_mu = D2{1}(mu(ind(:,i)));
209
210 p = ind(:,i);
211 obj.D2_lambda{1}(p,p) = D_lambda;
212 obj.D2_mu{1}(p,p) = D_mu;
213 end
214
215 for i = 1:m(1)
216 D_lambda = D2{2}(lambda(ind(i,:)));
217 D_mu = D2{2}(mu(ind(i,:)));
218
219 p = ind(i,:);
220 obj.D2_lambda{2}(p,p) = D_lambda;
221 obj.D2_mu{2}(p,p) = D_mu;
222 end
223
224 % Quadratures
225 obj.H = kron(H{1},H{2});
226 obj.Hi = inv(obj.H);
227 obj.H_w = H{2};
228 obj.H_e = H{2};
229 obj.H_s = H{1};
230 obj.H_n = H{1};
231 obj.H_1D = {H{1}, H{2}};
199 232
200 % Differentiation matrix D (without SAT) 233 % Differentiation matrix D (without SAT)
201 D2_lambda = obj.D2_lambda; 234 D2_lambda = obj.D2_lambda;
202 D2_mu = obj.D2_mu; 235 D2_mu = obj.D2_mu;
203 D1 = obj.D1; 236 D1 = obj.D1;
219 %=========================================%' 252 %=========================================%'
220 253
221 % Numerical traction operators for BC. 254 % Numerical traction operators for BC.
222 % Because d1 =/= e0^T*D1, the numerical tractions are different 255 % Because d1 =/= e0^T*D1, the numerical tractions are different
223 % at every boundary. 256 % at every boundary.
257 %
258 % Formula at boundary j: % tau^{j}_i = sum_k T^{j}_{ik} u_k
259 %
224 T_l = cell(dim,1); 260 T_l = cell(dim,1);
225 T_r = cell(dim,1); 261 T_r = cell(dim,1);
226 tau_l = cell(dim,1); 262 tau_l = cell(dim,1);
227 tau_r = cell(dim,1); 263 tau_r = cell(dim,1);
228 % tau^{j}_i = sum_k T^{j}_{ik} u_k 264
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; 265 D1 = obj.D1;
235 266
236 % Loop over boundaries 267 % Loop over boundaries
237 for j = 1:dim 268 for j = 1:dim
238 T_l{j} = cell(dim,dim); 269 T_l{j} = cell(dim,dim);
248 [~, n_l] = size(e_l{j}); 279 [~, n_l] = size(e_l{j});
249 [~, n_r] = size(e_r{j}); 280 [~, n_r] = size(e_r{j});
250 281
251 % Loop over components 282 % Loop over components
252 for i = 1:dim 283 for i = 1:dim
253 tau_l{j}{i} = sparse(n_l, dim*m_tot); 284 tau_l{j}{i} = sparse(dim*m_tot, n_l);
254 tau_r{j}{i} = sparse(n_r, dim*m_tot); 285 tau_r{j}{i} = sparse(dim*m_tot, n_r);
255 for k = 1:dim 286 for k = 1:dim
256 T_l{j}{i,k} = ... 287 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})... 288 (-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})... 289 -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}'; 290 -d(i,k)*MU_l*d1_l{j}')';
260 291
261 T_r{j}{i,k} = ... 292 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})... 293 (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})... 294 +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}'; 295 +d(i,k)*MU_r*d1_r{j}')';
265 296
266 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; 297 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}'; 298 tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,k}'*E{k}')';
268 end 299 end
269 300
270 end 301 end
271 end 302 end
272 303
273 % Transpose T and tau to match boundary operator convention 304 % Traction tensors, T_ij
274 for i = 1:dim 305 obj.T_w = T_l{1};
275 for j = 1:dim 306 obj.T_e = T_r{1};
276 tau_l{i}{j} = transpose(tau_l{i}{j}); 307 obj.T_s = T_l{2};
277 tau_r{i}{j} = transpose(tau_r{i}{j}); 308 obj.T_n = T_r{2};
278 for k = 1:dim 309
279 T_l{i}{j,k} = transpose(T_l{i}{j,k}); 310 % Restriction operators
280 T_r{i}{j,k} = transpose(T_r{i}{j,k}); 311 obj.e_w = e_w;
281 end 312 obj.e_e = e_e;
282 end 313 obj.e_s = e_s;
283 end 314 obj.e_n = e_n;
284 315
285 obj.T_l = T_l; 316 obj.e1_w = e1_w;
286 obj.T_r = T_r; 317 obj.e1_e = e1_e;
287 obj.tau_l = tau_l; 318 obj.e1_s = e1_s;
288 obj.tau_r = tau_r; 319 obj.e1_n = e1_n;
320
321 obj.e2_w = e2_w;
322 obj.e2_e = e2_e;
323 obj.e2_s = e2_s;
324 obj.e2_n = e2_n;
325
326 obj.e_scalar_w = e_scalar_w;
327 obj.e_scalar_e = e_scalar_e;
328 obj.e_scalar_s = e_scalar_s;
329 obj.e_scalar_n = e_scalar_n;
330
331 % First component of traction
332 obj.tau1_w = tau_l{1}{1};
333 obj.tau1_e = tau_r{1}{1};
334 obj.tau1_s = tau_l{2}{1};
335 obj.tau1_n = tau_r{2}{1};
336
337 % Second component of traction
338 obj.tau2_w = tau_l{1}{2};
339 obj.tau2_e = tau_r{1}{2};
340 obj.tau2_s = tau_l{2}{2};
341 obj.tau2_n = tau_r{2}{2};
342
343 % Traction vectors
344 obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')';
345 obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')';
346 obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')';
347 obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
289 348
290 % Kroneckered norms and coefficients 349 % Kroneckered norms and coefficients
291 I_dim = speye(dim);
292 obj.RHOi_kron = kron(obj.RHOi, I_dim); 350 obj.RHOi_kron = kron(obj.RHOi, I_dim);
293 obj.Hi_kron = kron(obj.Hi, I_dim); 351 obj.Hi_kron = kron(obj.Hi, I_dim);
294 352
295 % Misc. 353 % Misc.
296 obj.m = m; 354 obj.m = m;
298 obj.order = order; 356 obj.order = order;
299 obj.grid = g; 357 obj.grid = g;
300 obj.dim = dim; 358 obj.dim = dim;
301 359
302 % B, used for adjoint optimization 360 % B, used for adjoint optimization
303 B = cell(dim, 1); 361 B = [];
304 for i = 1:dim 362 if optFlag
305 B{i} = cell(m_tot, 1); 363 B = cell(dim, 1);
306 end 364 for i = 1:dim
307 365 B{i} = cell(m_tot, 1);
308 for i = 1:dim 366 end
309 for j = 1:m_tot 367
310 B{i}{j} = sparse(m_tot, m_tot); 368 B0 = sparse(m_tot, m_tot);
311 end 369 for i = 1:dim
312 end 370 for j = 1:m_tot
313 371 B{i}{j} = B0;
314 ind = grid.funcToMatrix(g, 1:m_tot); 372 end
315 373 end
316 % Direction 1 374
317 for k = 1:m(1) 375 ind = grid.funcToMatrix(g, 1:m_tot);
318 c = sparse(m(1),1); 376
319 c(k) = 1; 377 % Direction 1
320 [~, B_1D] = ops{1}.D2(c); 378 for k = 1:m(1)
321 for l = 1:m(2) 379 c = sparse(m(1),1);
322 p = ind(:,l); 380 c(k) = 1;
323 B{1}{(k-1)*m(2) + l}(p, p) = B_1D; 381 [~, B_1D] = ops{1}.D2(c);
324 end 382 for l = 1:m(2)
325 end 383 p = ind(:,l);
326 384 B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
327 % Direction 2 385 end
328 for k = 1:m(2) 386 end
329 c = sparse(m(2),1); 387
330 c(k) = 1; 388 % Direction 2
331 [~, B_1D] = ops{2}.D2(c); 389 for k = 1:m(2)
332 for l = 1:m(1) 390 c = sparse(m(2),1);
333 p = ind(l,:); 391 c(k) = 1;
334 B{2}{(l-1)*m(2) + k}(p, p) = B_1D; 392 [~, B_1D] = ops{2}.D2(c);
335 end 393 for l = 1:m(1)
336 end 394 p = ind(l,:);
337 395 B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
396 end
397 end
398 end
338 obj.B = B; 399 obj.B = B;
400
339 401
340 end 402 end
341 403
342 404
343 % Closure functions return the operators applied to the own domain to close the boundary 405 % Closure functions return the operators applied to the own domain to close the boundary
344 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. 406 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
345 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. 407 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
346 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition 408 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
347 % on the first component. 409 % on the first component. Can also be e.g.
410 % {'normal', 'd'} or {'tangential', 't'} for conditions on
411 % tangential/normal component.
348 % data is a function returning the data that should be applied at the boundary. 412 % data is a function returning the data that should be applied at the boundary.
349 % neighbour_scheme is an instance of Scheme that should be interfaced to. 413 % neighbour_scheme is an instance of Scheme that should be interfaced to.
350 % neighbour_boundary is a string specifying which boundary to interface to. 414 % neighbour_boundary is a string specifying which boundary to interface to.
351 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) 415 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
352 default_arg('tuning', 1.2); 416 default_arg('tuning', 1.2);
353 417
354 assert( iscell(bc), 'The BC type must be a 2x1 cell array' ); 418 assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
355 comp = bc{1}; 419 comp = bc{1};
356 type = bc{2}; 420 type = bc{2};
357 421 if ischar(comp)
358 % j is the coordinate direction of the boundary 422 comp = obj.getComponent(comp, boundary);
359 j = obj.get_boundary_number(boundary); 423 end
360 [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); 424
361 425 e = obj.getBoundaryOperatorForScalarField('e', boundary);
426 tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
427 T = obj.getBoundaryTractionOperator(boundary);
428 alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary);
429 H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
362 430
363 E = obj.E; 431 E = obj.E;
364 Hi = obj.Hi; 432 Hi = obj.Hi;
365 LAMBDA = obj.LAMBDA; 433 LAMBDA = obj.LAMBDA;
366 MU = obj.MU; 434 MU = obj.MU;
368 436
369 dim = obj.dim; 437 dim = obj.dim;
370 m_tot = obj.grid.N(); 438 m_tot = obj.grid.N();
371 439
372 % Preallocate 440 % Preallocate
441 [~, col] = size(tau);
373 closure = sparse(dim*m_tot, dim*m_tot); 442 closure = sparse(dim*m_tot, dim*m_tot);
374 penalty = sparse(dim*m_tot, m_tot/obj.m(j)); 443 penalty = sparse(dim*m_tot, col);
375 444
376 k = comp; 445 k = comp;
377 switch type 446 switch type
378 447
379 % Dirichlet boundary condition 448 % Dirichlet boundary condition
380 case {'D','d','dirichlet','Dirichlet'} 449 case {'D','d','dirichlet','Dirichlet'}
381
382 alpha = obj.getBoundaryOperator('alpha', boundary);
383 450
384 % Loop over components that Dirichlet penalties end up on 451 % Loop over components that Dirichlet penalties end up on
385 for i = 1:dim 452 for i = 1:dim
386 C = transpose(T{k,i}); 453 C = transpose(T{k,i});
387 A = -tuning*e*transpose(alpha{i,k}); 454 A = -tuning*e*transpose(alpha{i,k});
390 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; 457 penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
391 end 458 end
392 459
393 % Free boundary condition 460 % Free boundary condition
394 case {'F','f','Free','free','traction','Traction','t','T'} 461 case {'F','f','Free','free','traction','Traction','t','T'}
395 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}'; 462 closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau';
396 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; 463 penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
397 464
398 % Unknown boundary condition 465 % Unknown boundary condition
399 otherwise 466 otherwise
400 error('No such boundary condition: type = %s',type); 467 error('No such boundary condition: type = %s',type);
427 % u denotes the solution in the own domain 494 % u denotes the solution in the own domain
428 % v denotes the solution in the neighbour domain 495 % v denotes the solution in the neighbour domain
429 % Operators without subscripts are from the own domain. 496 % Operators without subscripts are from the own domain.
430 497
431 % Get boundary operators 498 % Get boundary operators
432 e = obj.getBoundaryOperator('e_tot', boundary); 499 e = obj.getBoundaryOperator('e', boundary);
433 tau = obj.getBoundaryOperator('tau_tot', boundary); 500 tau = obj.getBoundaryOperator('tau', boundary);
434 501
435 e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary); 502 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
436 tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary); 503 tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary);
437 504
438 H_gamma = obj.getBoundaryQuadrature(boundary); 505 H_gamma = obj.getBoundaryQuadrature(boundary);
439 506
440 % Operators and quantities that correspond to the own domain only 507 % Operators and quantities that correspond to the own domain only
441 Hi = obj.Hi_kron; 508 Hi = obj.Hi_kron;
442 RHOi = obj.RHOi_kron; 509 RHOi = obj.RHOi_kron;
443 510
444 % Penalty strength operators 511 % Penalty strength operators
445 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary); 512 alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary);
446 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary); 513 alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary);
447 514
448 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e'); 515 closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
449 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v'); 516 penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
450 517
451 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau'; 518 closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
458 525
459 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) 526 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
460 error('Non-conforming interfaces not implemented yet.'); 527 error('Non-conforming interfaces not implemented yet.');
461 end 528 end
462 529
463 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 530 % Returns the component number that is the tangential/normal component
464 function [j, nj] = get_boundary_number(obj, boundary) 531 % at the specified boundary
465 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 532 function comp = getComponent(obj, comp_str, boundary)
533 assertIsMember(comp_str, {'normal', 'tangential'});
534 assertIsMember(boundary, {'w', 'e', 's', 'n'});
466 535
467 switch boundary 536 switch boundary
468 case {'w', 'e'} 537 case {'w', 'e'}
469 j = 1; 538 switch comp_str
470 case {'s', 'n'} 539 case 'normal'
471 j = 2; 540 comp = 1;
472 end 541 case 'tangential'
473 542 comp = 2;
474 switch boundary 543 end
475 case {'w', 's'} 544 case {'s', 'n'}
476 nj = -1; 545 switch comp_str
477 case {'e', 'n'} 546 case 'normal'
478 nj = 1; 547 comp = 2;
548 case 'tangential'
549 comp = 1;
550 end
479 end 551 end
480 end 552 end
481 553
482 % Returns the boundary operator op for the boundary specified by the string boundary. 554 % Returns the boundary operator op for the boundary specified by the string boundary.
483 % op -- string 555 % op -- string
484 % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() 556 function o = getBoundaryOperator(obj, op, boundary)
485 function [varargout] = getBoundaryOperator(obj, op, boundary)
486 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 557 assertIsMember(boundary, {'w', 'e', 's', 'n'})
487 assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'}) 558 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha', 'alpha1', 'alpha2'})
488
489 switch boundary
490 case {'w', 'e'}
491 j = 1;
492 case {'s', 'n'}
493 j = 2;
494 end
495 559
496 switch op 560 switch op
561
562 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
563 o = obj.([op, '_', boundary]);
564
565 % Yields vector-valued penalty strength given displacement BC on all components
566 case 'alpha'
567 e = obj.getBoundaryOperator('e', boundary);
568 e_scalar = obj.getBoundaryOperatorForScalarField('e', boundary);
569 alpha_scalar = obj.getBoundaryOperatorForScalarField('alpha', boundary);
570 E = obj.E;
571 [m, n] = size(alpha_scalar{1,1});
572 alpha = sparse(m*obj.dim, n*obj.dim);
573 for i = 1:obj.dim
574 for l = 1:obj.dim
575 alpha = alpha + (e'*E{i}*e_scalar*alpha_scalar{i,l}'*E{l}')';
576 end
577 end
578 o = alpha;
579
580 % Yields penalty strength for component 1 given displacement BC on all components
581 case 'alpha1'
582 alpha = obj.getBoundaryOperator('alpha', boundary);
583 e = obj.getBoundaryOperator('e', boundary);
584 e1 = obj.getBoundaryOperator('e1', boundary);
585
586 alpha1 = (e1'*e*alpha')';
587 o = alpha1;
588
589 % Yields penalty strength for component 2 given displacement BC on all components
590 case 'alpha2'
591 alpha = obj.getBoundaryOperator('alpha', boundary);
592 e = obj.getBoundaryOperator('e', boundary);
593 e2 = obj.getBoundaryOperator('e2', boundary);
594
595 alpha2 = (e2'*e*alpha')';
596 o = alpha2;
597 end
598
599 end
600
601 % Returns the boundary operator op for the boundary specified by the string boundary.
602 % op -- string
603 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
604 assertIsMember(boundary, {'w', 'e', 's', 'n'})
605 assertIsMember(op, {'e', 'alpha'})
606
607 switch op
608
497 case 'e' 609 case 'e'
498 switch boundary 610 o = obj.(['e_scalar', '_', boundary]);
499 case {'w', 's'}
500 o = obj.e_l{j};
501 case {'e', 'n'}
502 o = obj.e_r{j};
503 end
504
505 case 'e_tot'
506 e = obj.getBoundaryOperator('e', boundary);
507 I_dim = speye(obj.dim, obj.dim);
508 o = kron(e, I_dim);
509
510 case 'd'
511 switch boundary
512 case {'w', 's'}
513 o = obj.d1_l{j};
514 case {'e', 'n'}
515 o = obj.d1_r{j};
516 end
517
518 case 'T'
519 switch boundary
520 case {'w', 's'}
521 o = obj.T_l{j};
522 case {'e', 'n'}
523 o = obj.T_r{j};
524 end
525
526 case 'tau'
527 switch boundary
528 case {'w', 's'}
529 o = obj.tau_l{j};
530 case {'e', 'n'}
531 o = obj.tau_r{j};
532 end
533
534 case 'tau_tot'
535 [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
536
537 I_dim = speye(obj.dim, obj.dim);
538 e_tot = kron(e, I_dim);
539 E = obj.E;
540 tau_tot = (e_tot'*E{1}*e*tau{1}')';
541 for i = 2:obj.dim
542 tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')';
543 end
544 o = tau_tot;
545
546 case 'H'
547 o = obj.H_boundary{j};
548 611
549 case 'alpha' 612 case 'alpha'
550 % alpha = alpha(i,j) is the penalty strength for displacement BC. 613
551 e = obj.getBoundaryOperator('e', boundary); 614 % alpha{i,j} is the penalty strength on component i due to
615 % displacement BC for component j.
616 e = obj.getBoundaryOperatorForScalarField('e', boundary);
552 617
553 LAMBDA = obj.LAMBDA; 618 LAMBDA = obj.LAMBDA;
554 MU = obj.MU; 619 MU = obj.MU;
555
556 dim = obj.dim; 620 dim = obj.dim;
557 theta_R = obj.theta_R{j}; 621
558 theta_H = obj.theta_H{j}; 622 switch boundary
559 theta_M = obj.theta_M{j}; 623 case {'w', 'e'}
624 k = 1;
625 case {'s', 'n'}
626 k = 2;
627 end
628
629 theta_R = obj.theta_R{k};
630 theta_H = obj.theta_H{k};
631 theta_M = obj.theta_M{k};
560 632
561 a_lambda = dim/theta_H + 1/theta_R; 633 a_lambda = dim/theta_H + 1/theta_R;
562 a_mu_i = 2/theta_M; 634 a_mu_i = 2/theta_M;
563 a_mu_ij = 2/theta_H + 1/theta_R; 635 a_mu_ij = 2/theta_H + 1/theta_R;
564 636
565 d = @kroneckerDelta; % Kronecker delta 637 d = @kroneckerDelta; % Kronecker delta
566 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta 638 db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
567 alpha = cell(obj.dim, obj.dim);
568 639
569 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ... 640 alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ...
570 + d(i,j)* a_mu_i*MU ... 641 + d(i,j)* a_mu_i*MU ...
571 + db(i,j)*a_mu_ij*MU; 642 + db(i,j)*a_mu_ij*MU;
643
644 alpha = cell(obj.dim, obj.dim);
572 for i = 1:obj.dim 645 for i = 1:obj.dim
573 for l = 1:obj.dim 646 for j = 1:obj.dim
574 alpha{i,l} = d(i,l)*alpha_func(i,j)*e; 647 alpha{i,j} = d(i,j)*alpha_func(i,k)*e;
575 end 648 end
576 end 649 end
577
578 o = alpha; 650 o = alpha;
579 651 end
580 case 'alpha_tot' 652
581 % alpha = alpha(i,j) is the penalty strength for displacement BC. 653 end
582 [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary); 654
583 E = obj.E; 655 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
584 [m, n] = size(alpha{1,1}); 656 % Formula: tau_i = T_ij u_j
585 alpha_tot = sparse(m*obj.dim, n*obj.dim); 657 % op -- string
586 for i = 1:obj.dim 658 function T = getBoundaryTractionOperator(obj, boundary)
587 for l = 1:obj.dim 659 assertIsMember(boundary, {'w', 'e', 's', 'n'})
588 alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; 660
589 end 661 T = obj.(['T', '_', boundary]);
590 end
591 o = alpha_tot;
592 end
593
594 end 662 end
595 663
596 % Returns square boundary quadrature matrix, of dimension 664 % Returns square boundary quadrature matrix, of dimension
597 % corresponding to the number of boundary points 665 % corresponding to the number of boundary unknowns
598 % 666 %
599 % boundary -- string 667 % boundary -- string
600 function H = getBoundaryQuadrature(obj, boundary) 668 function H = getBoundaryQuadrature(obj, boundary)
601 assertIsMember(boundary, {'w', 'e', 's', 'n'}) 669 assertIsMember(boundary, {'w', 'e', 's', 'n'})
602 670
603 switch boundary 671 H = obj.getBoundaryQuadratureForScalarField(boundary);
604 case {'w','e'}
605 j = 1;
606 case {'s','n'}
607 j = 2;
608 end
609 H = obj.H_boundary{j};
610 I_dim = speye(obj.dim, obj.dim); 672 I_dim = speye(obj.dim, obj.dim);
611 H = kron(H, I_dim); 673 H = kron(H, I_dim);
612 end 674 end
613 675
676 % Returns square boundary quadrature matrix, of dimension
677 % corresponding to the number of boundary grid points
678 %
679 % boundary -- string
680 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
681 assertIsMember(boundary, {'w', 'e', 's', 'n'})
682
683 H_b = obj.(['H_', boundary]);
684 end
685
614 function N = size(obj) 686 function N = size(obj)
615 N = obj.dim*prod(obj.m); 687 N = obj.dim*prod(obj.m);
616 end 688 end
617 end 689 end
618 end 690 end