Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropicUpwind.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 8ff3a95ad7cc
children
comparison
equal deleted inserted replaced
1330:855871e0b852 1331:60c875c18de3
1 classdef Elastic2dCurvilinearAnisotropicUpwind < 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 Dx, Dy % Physical derivatives
27 sigma % Cell matrix of physical stress operators
28 n_w, n_e, n_s, n_n % Physical normals
29
30 % Boundary operators in cell format, used for BC
31 T_w, T_e, T_s, T_n
32
33 % Traction operators
34 tau_w, tau_e, tau_s, tau_n % Return vector field
35 tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field
36 tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field
37
38 % Inner products
39 H
40
41 % Boundary inner products (for scalar field)
42 H_w, H_e, H_s, H_n
43
44 % Surface Jacobian vectors
45 s_w, s_e, s_s, s_n
46
47 % Boundary restriction operators
48 e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary
49 e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary
50 e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary
51 e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
52 en_w, en_e, en_s, en_n % Act on vector field, return normal component
53
54 % E{i}^T picks out component i
55 E
56
57 % Elastic2dVariableAnisotropic object for reference domain
58 refObj
59 end
60
61 methods
62
63 % The coefficients can either be function handles or grid functions
64 % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
65 function obj = Elastic2dCurvilinearAnisotropicUpwind(g, order, rho, C, opSet, optFlag)
66 default_arg('rho', @(x,y) 0*x+1);
67 default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind});
68 default_arg('optFlag', false);
69 dim = 2;
70
71 C_default = cell(dim,dim,dim,dim);
72 for i = 1:dim
73 for j = 1:dim
74 for k = 1:dim
75 for l = 1:dim
76 C_default{i,j,k,l} = @(x,y) 0*x;
77 end
78 end
79 end
80 end
81 default_arg('C', C_default);
82
83 assert(isa(g, 'grid.Curvilinear'));
84
85 if isa(rho, 'function_handle')
86 rho = grid.evalOn(g, rho);
87 end
88
89 C_mat = cell(dim,dim,dim,dim);
90 for i = 1:dim
91 for j = 1:dim
92 for k = 1:dim
93 for l = 1:dim
94 if isa(C{i,j,k,l}, 'function_handle')
95 C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l});
96 end
97 C_mat{i,j,k,l} = spdiag(C{i,j,k,l});
98 end
99 end
100 end
101 end
102 obj.C = C_mat;
103
104 m = g.size();
105 m_tot = g.N();
106
107 % 1D operators
108 opSetMetric = {@sbp.D2VariableCompatible, @sbp.D2VariableCompatible};
109 orderMetric = ceil(order/2)*2;
110 m_u = m(1);
111 m_v = m(2);
112 ops_u = opSetMetric{1}(m_u, {0, 1}, orderMetric);
113 ops_v = opSetMetric{2}(m_v, {0, 1}, orderMetric);
114
115 h_u = ops_u.h;
116 h_v = ops_v.h;
117
118 I_u = speye(m_u);
119 I_v = speye(m_v);
120
121 D1_u = ops_u.D1;
122 H_u = ops_u.H;
123 Hi_u = ops_u.HI;
124 e_l_u = ops_u.e_l;
125 e_r_u = ops_u.e_r;
126 d1_l_u = ops_u.d1_l;
127 d1_r_u = ops_u.d1_r;
128
129 D1_v = ops_v.D1;
130 H_v = ops_v.H;
131 Hi_v = ops_v.HI;
132 e_l_v = ops_v.e_l;
133 e_r_v = ops_v.e_r;
134 d1_l_v = ops_v.d1_l;
135 d1_r_v = ops_v.d1_r;
136
137
138 % Logical operators
139 Du = kr(D1_u,I_v);
140 Dv = kr(I_u,D1_v);
141
142 e_w = kr(e_l_u,I_v);
143 e_e = kr(e_r_u,I_v);
144 e_s = kr(I_u,e_l_v);
145 e_n = kr(I_u,e_r_v);
146
147 % Metric coefficients
148 coords = g.points();
149 x = coords(:,1);
150 y = coords(:,2);
151
152 x_u = Du*x;
153 x_v = Dv*x;
154 y_u = Du*y;
155 y_v = Dv*y;
156
157 J = x_u.*y_v - x_v.*y_u;
158
159 K = cell(dim, dim);
160 K{1,1} = y_v./J;
161 K{1,2} = -y_u./J;
162 K{2,1} = -x_v./J;
163 K{2,2} = x_u./J;
164
165 % Physical derivatives
166 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
167 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
168
169 % Wrap around Aniosotropic Cartesian
170 rho_tilde = J.*rho;
171
172 PHI = cell(dim,dim,dim,dim);
173 for i = 1:dim
174 for j = 1:dim
175 for k = 1:dim
176 for l = 1:dim
177 PHI{i,j,k,l} = 0*C{i,j,k,l};
178 for m = 1:dim
179 for n = 1:dim
180 PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,i}.*C{m,j,n,l}.*K{n,k};
181 end
182 end
183 end
184 end
185 end
186 end
187
188 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
189 refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, opSet);
190
191 %---- Set object properties ------
192 obj.RHO = spdiag(rho);
193
194 % Volume quadrature
195 obj.J = spdiag(J);
196 obj.Ji = spdiag(1./J);
197 obj.H = obj.J*refObj.H;
198
199 % Boundary quadratures
200 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
201 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
202 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
203 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
204 obj.s_w = s_w;
205 obj.s_e = s_e;
206 obj.s_s = s_s;
207 obj.s_n = s_n;
208
209 obj.H_w = H_v*spdiag(s_w);
210 obj.H_e = H_v*spdiag(s_e);
211 obj.H_s = H_u*spdiag(s_s);
212 obj.H_n = H_u*spdiag(s_n);
213
214 % Restriction operators
215 obj.e_w = refObj.e_w;
216 obj.e_e = refObj.e_e;
217 obj.e_s = refObj.e_s;
218 obj.e_n = refObj.e_n;
219
220 % Adapt things from reference object
221 obj.D = refObj.D;
222 obj.E = refObj.E;
223
224 obj.e1_w = refObj.e1_w;
225 obj.e1_e = refObj.e1_e;
226 obj.e1_s = refObj.e1_s;
227 obj.e1_n = refObj.e1_n;
228
229 obj.e2_w = refObj.e2_w;
230 obj.e2_e = refObj.e2_e;
231 obj.e2_s = refObj.e2_s;
232 obj.e2_n = refObj.e2_n;
233
234 obj.e_scalar_w = refObj.e_scalar_w;
235 obj.e_scalar_e = refObj.e_scalar_e;
236 obj.e_scalar_s = refObj.e_scalar_s;
237 obj.e_scalar_n = refObj.e_scalar_n;
238
239 e1_w = obj.e1_w;
240 e1_e = obj.e1_e;
241 e1_s = obj.e1_s;
242 e1_n = obj.e1_n;
243
244 e2_w = obj.e2_w;
245 e2_e = obj.e2_e;
246 e2_s = obj.e2_s;
247 e2_n = obj.e2_n;
248
249 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
250 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
251 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
252 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
253
254 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
255 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
256 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
257 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
258
259 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
260 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
261 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
262 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
263
264 % Physical normals
265 e_w = obj.e_scalar_w;
266 e_e = obj.e_scalar_e;
267 e_s = obj.e_scalar_s;
268 e_n = obj.e_scalar_n;
269
270 e_w_vec = obj.e_w;
271 e_e_vec = obj.e_e;
272 e_s_vec = obj.e_s;
273 e_n_vec = obj.e_n;
274
275 nu_w = [-1,0];
276 nu_e = [1,0];
277 nu_s = [0,-1];
278 nu_n = [0,1];
279
280 obj.n_w = cell(2,1);
281 obj.n_e = cell(2,1);
282 obj.n_s = cell(2,1);
283 obj.n_n = cell(2,1);
284
285 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
286 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
287 obj.n_w{1} = spdiag(n_w_1);
288 obj.n_w{2} = spdiag(n_w_2);
289
290 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
291 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
292 obj.n_e{1} = spdiag(n_e_1);
293 obj.n_e{2} = spdiag(n_e_2);
294
295 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
296 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
297 obj.n_s{1} = spdiag(n_s_1);
298 obj.n_s{2} = spdiag(n_s_2);
299
300 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
301 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
302 obj.n_n{1} = spdiag(n_n_1);
303 obj.n_n{2} = spdiag(n_n_2);
304
305 % Operators that extract the normal component
306 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
307 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
308 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
309 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
310
311 % Stress operators
312 sigma = cell(dim, dim);
313 D1 = {obj.Dx, obj.Dy};
314 E = obj.E;
315 N = length(obj.RHO);
316 for i = 1:dim
317 for j = 1:dim
318 sigma{i,j} = sparse(N,2*N);
319 for k = 1:dim
320 for l = 1:dim
321 sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}';
322 end
323 end
324 end
325 end
326 obj.sigma = sigma;
327
328 % Misc.
329 obj.refObj = refObj;
330 obj.m = refObj.m;
331 obj.h = refObj.h;
332 obj.order = order;
333 obj.grid = g;
334 obj.dim = dim;
335
336 end
337
338
339 % Closure functions return the operators applied to the own domain to close the boundary
340 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
341 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
342 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
343 % on the first component. Can also be e.g.
344 % {'normal', 'd'} or {'tangential', 't'} for conditions on
345 % tangential/normal component.
346 % data is a function returning the data that should be applied at the boundary.
347 % neighbour_scheme is an instance of Scheme that should be interfaced to.
348 % neighbour_boundary is a string specifying which boundary to interface to.
349
350 % For displacement bc:
351 % bc = {comp, 'd', dComps},
352 % where
353 % dComps = vector of components with displacement BC. Default: 1:dim.
354 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
355 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
356 default_arg('tuning', 1.0);
357 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
358
359 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
360
361 type = bc{2};
362
363 switch type
364 case {'F','f','Free','free','traction','Traction','t','T'}
365 s = obj.(['s_' boundary]);
366 s = spdiag(s);
367 penalty = penalty*s;
368 end
369 end
370
371 % type Struct that specifies the interface coupling.
372 % Fields:
373 % -- tuning: penalty strength, defaults to 1.0
374 % -- interpolation: type of interpolation, default 'none'
375 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
376
377 defaultType.tuning = 1.0;
378 defaultType.interpolation = 'none';
379 default_struct('type', defaultType);
380
381 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
382 end
383
384 % Returns h11 for the boundary specified by the string boundary.
385 % op -- string
386 function h11 = getBorrowing(obj, boundary)
387 assertIsMember(boundary, {'w', 'e', 's', 'n'})
388
389 switch boundary
390 case {'w','e'}
391 h11 = obj.refObj.h11{1};
392 case {'s', 'n'}
393 h11 = obj.refObj.h11{2};
394 end
395 end
396
397 % Returns the outward unit normal vector for the boundary specified by the string boundary.
398 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
399 function n = getNormal(obj, boundary)
400 assertIsMember(boundary, {'w', 'e', 's', 'n'})
401
402 n = obj.(['n_' boundary]);
403 end
404
405 % Returns the boundary operator op for the boundary specified by the string boundary.
406 % op -- string
407 function o = getBoundaryOperator(obj, op, boundary)
408 assertIsMember(boundary, {'w', 'e', 's', 'n'})
409 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
410
411 o = obj.([op, '_', boundary]);
412
413 end
414
415 % Returns the boundary operator op for the boundary specified by the string boundary.
416 % op -- string
417 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
418 assertIsMember(boundary, {'w', 'e', 's', 'n'})
419 assertIsMember(op, {'e'})
420
421 switch op
422
423 case 'e'
424 o = obj.(['e_scalar', '_', boundary]);
425 end
426
427 end
428
429 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
430 % Formula: tau_i = T_ij u_j
431 % op -- string
432 function T = getBoundaryTractionOperator(obj, boundary)
433 assertIsMember(boundary, {'w', 'e', 's', 'n'})
434
435 T = obj.(['T', '_', boundary]);
436 end
437
438 % Returns square boundary quadrature matrix, of dimension
439 % corresponding to the number of boundary unknowns
440 %
441 % boundary -- string
442 function H = getBoundaryQuadrature(obj, boundary)
443 assertIsMember(boundary, {'w', 'e', 's', 'n'})
444
445 H = obj.getBoundaryQuadratureForScalarField(boundary);
446 I_dim = speye(obj.dim, obj.dim);
447 H = kron(H, I_dim);
448 end
449
450 % Returns square boundary quadrature matrix, of dimension
451 % corresponding to the number of boundary grid points
452 %
453 % boundary -- string
454 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
455 assertIsMember(boundary, {'w', 'e', 's', 'n'})
456
457 H_b = obj.(['H_', boundary]);
458 end
459
460 function N = size(obj)
461 N = obj.dim*prod(obj.m);
462 end
463 end
464 end