Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropicUpwind.m @ 1227:02dfe3a56742 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Upwind ElasticAnisotropic schemes. Seem to work really well!
author Martin Almquist <malmquist@stanford.edu>
date Sat, 16 Nov 2019 14:26:06 -0800
parents
children 9a7cac4202d0
comparison
equal deleted inserted replaced
1214:e1d4cb8b5309 1227:02dfe3a56742
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.D2VariableCompatible, @sbp.D2VariableCompatible});
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 + 1;
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 m_u = m(1);
109 m_v = m(2);
110 ops_u = opSet{1}(m_u, {0, 1}, order);
111 ops_v = opSet{2}(m_v, {0, 1}, order);
112
113 h_u = ops_u.h;
114 h_v = ops_v.h;
115
116 I_u = speye(m_u);
117 I_v = speye(m_v);
118
119 D1_u = ops_u.D1;
120 H_u = ops_u.H;
121 Hi_u = ops_u.HI;
122 e_l_u = ops_u.e_l;
123 e_r_u = ops_u.e_r;
124 d1_l_u = ops_u.d1_l;
125 d1_r_u = ops_u.d1_r;
126
127 D1_v = ops_v.D1;
128 H_v = ops_v.H;
129 Hi_v = ops_v.HI;
130 e_l_v = ops_v.e_l;
131 e_r_v = ops_v.e_r;
132 d1_l_v = ops_v.d1_l;
133 d1_r_v = ops_v.d1_r;
134
135
136 % Logical operators
137 Du = kr(D1_u,I_v);
138 Dv = kr(I_u,D1_v);
139
140 e_w = kr(e_l_u,I_v);
141 e_e = kr(e_r_u,I_v);
142 e_s = kr(I_u,e_l_v);
143 e_n = kr(I_u,e_r_v);
144
145 % Metric coefficients
146 coords = g.points();
147 x = coords(:,1);
148 y = coords(:,2);
149
150 x_u = Du*x;
151 x_v = Dv*x;
152 y_u = Du*y;
153 y_v = Dv*y;
154
155 J = x_u.*y_v - x_v.*y_u;
156
157 K = cell(dim, dim);
158 K{1,1} = y_v./J;
159 K{1,2} = -y_u./J;
160 K{2,1} = -x_v./J;
161 K{2,2} = x_u./J;
162
163 % Physical derivatives
164 obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv;
165 obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv;
166
167 % Wrap around Aniosotropic Cartesian
168 rho_tilde = J.*rho;
169
170 PHI = cell(dim,dim,dim,dim);
171 for i = 1:dim
172 for j = 1:dim
173 for k = 1:dim
174 for l = 1:dim
175 PHI{i,j,k,l} = 0*C{i,j,k,l};
176 for m = 1:dim
177 for n = 1:dim
178 PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,i}.*C{m,j,n,l}.*K{n,k};
179 end
180 end
181 end
182 end
183 end
184 end
185
186 gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1});
187 refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, []);
188
189 %---- Set object properties ------
190 obj.RHO = spdiag(rho);
191
192 % Volume quadrature
193 obj.J = spdiag(J);
194 obj.Ji = spdiag(1./J);
195 obj.H = obj.J*kr(H_u,H_v);
196
197 % Boundary quadratures
198 s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2);
199 s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2);
200 s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2);
201 s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2);
202 obj.s_w = s_w;
203 obj.s_e = s_e;
204 obj.s_s = s_s;
205 obj.s_n = s_n;
206
207 obj.H_w = H_v*spdiag(s_w);
208 obj.H_e = H_v*spdiag(s_e);
209 obj.H_s = H_u*spdiag(s_s);
210 obj.H_n = H_u*spdiag(s_n);
211
212 % Restriction operators
213 obj.e_w = refObj.e_w;
214 obj.e_e = refObj.e_e;
215 obj.e_s = refObj.e_s;
216 obj.e_n = refObj.e_n;
217
218 % Adapt things from reference object
219 obj.D = refObj.D;
220 obj.E = refObj.E;
221
222 obj.e1_w = refObj.e1_w;
223 obj.e1_e = refObj.e1_e;
224 obj.e1_s = refObj.e1_s;
225 obj.e1_n = refObj.e1_n;
226
227 obj.e2_w = refObj.e2_w;
228 obj.e2_e = refObj.e2_e;
229 obj.e2_s = refObj.e2_s;
230 obj.e2_n = refObj.e2_n;
231
232 obj.e_scalar_w = refObj.e_scalar_w;
233 obj.e_scalar_e = refObj.e_scalar_e;
234 obj.e_scalar_s = refObj.e_scalar_s;
235 obj.e_scalar_n = refObj.e_scalar_n;
236
237 e1_w = obj.e1_w;
238 e1_e = obj.e1_e;
239 e1_s = obj.e1_s;
240 e1_n = obj.e1_n;
241
242 e2_w = obj.e2_w;
243 e2_e = obj.e2_e;
244 e2_s = obj.e2_s;
245 e2_n = obj.e2_n;
246
247 obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
248 obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
249 obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')';
250 obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')';
251
252 obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')';
253 obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')';
254 obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')';
255 obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')';
256
257 obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')';
258 obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')';
259 obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
260 obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
261
262 % Physical normals
263 e_w = obj.e_scalar_w;
264 e_e = obj.e_scalar_e;
265 e_s = obj.e_scalar_s;
266 e_n = obj.e_scalar_n;
267
268 e_w_vec = obj.e_w;
269 e_e_vec = obj.e_e;
270 e_s_vec = obj.e_s;
271 e_n_vec = obj.e_n;
272
273 nu_w = [-1,0];
274 nu_e = [1,0];
275 nu_s = [0,-1];
276 nu_n = [0,1];
277
278 obj.n_w = cell(2,1);
279 obj.n_e = cell(2,1);
280 obj.n_s = cell(2,1);
281 obj.n_n = cell(2,1);
282
283 n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
284 n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
285 obj.n_w{1} = spdiag(n_w_1);
286 obj.n_w{2} = spdiag(n_w_2);
287
288 n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
289 n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
290 obj.n_e{1} = spdiag(n_e_1);
291 obj.n_e{2} = spdiag(n_e_2);
292
293 n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
294 n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
295 obj.n_s{1} = spdiag(n_s_1);
296 obj.n_s{2} = spdiag(n_s_2);
297
298 n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
299 n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
300 obj.n_n{1} = spdiag(n_n_1);
301 obj.n_n{2} = spdiag(n_n_2);
302
303 % Operators that extract the normal component
304 obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
305 obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
306 obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
307 obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
308
309 % Stress operators
310 sigma = cell(dim, dim);
311 D1 = {obj.Dx, obj.Dy};
312 E = obj.E;
313 N = length(obj.RHO);
314 for i = 1:dim
315 for j = 1:dim
316 sigma{i,j} = sparse(N,2*N);
317 for k = 1:dim
318 for l = 1:dim
319 sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}';
320 end
321 end
322 end
323 end
324 obj.sigma = sigma;
325
326 % Misc.
327 obj.refObj = refObj;
328 obj.m = refObj.m;
329 obj.h = refObj.h;
330 obj.order = order;
331 obj.grid = g;
332 obj.dim = dim;
333
334 end
335
336
337 % Closure functions return the operators applied to the own domain to close the boundary
338 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
339 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
340 % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
341 % on the first component. Can also be e.g.
342 % {'normal', 'd'} or {'tangential', 't'} for conditions on
343 % tangential/normal component.
344 % data is a function returning the data that should be applied at the boundary.
345 % neighbour_scheme is an instance of Scheme that should be interfaced to.
346 % neighbour_boundary is a string specifying which boundary to interface to.
347
348 % For displacement bc:
349 % bc = {comp, 'd', dComps},
350 % where
351 % dComps = vector of components with displacement BC. Default: 1:dim.
352 % In this way, we can specify one BC at a time even though the SATs depend on all BC.
353 function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
354 default_arg('tuning', 1.0);
355 assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
356
357 [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
358
359 type = bc{2};
360
361 switch type
362 case {'F','f','Free','free','traction','Traction','t','T'}
363 s = obj.(['s_' boundary]);
364 s = spdiag(s);
365 penalty = penalty*s;
366 end
367 end
368
369 % type Struct that specifies the interface coupling.
370 % Fields:
371 % -- tuning: penalty strength, defaults to 1.0
372 % -- interpolation: type of interpolation, default 'none'
373 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
374
375 defaultType.tuning = 1.0;
376 defaultType.interpolation = 'none';
377 default_struct('type', defaultType);
378
379 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
380 end
381
382 % Returns h11 for the boundary specified by the string boundary.
383 % op -- string
384 function h11 = getBorrowing(obj, boundary)
385 assertIsMember(boundary, {'w', 'e', 's', 'n'})
386
387 switch boundary
388 case {'w','e'}
389 h11 = obj.refObj.h11{1};
390 case {'s', 'n'}
391 h11 = obj.refObj.h11{2};
392 end
393 end
394
395 % Returns the outward unit normal vector for the boundary specified by the string boundary.
396 % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
397 function n = getNormal(obj, boundary)
398 assertIsMember(boundary, {'w', 'e', 's', 'n'})
399
400 n = obj.(['n_' boundary]);
401 end
402
403 % Returns the boundary operator op for the boundary specified by the string boundary.
404 % op -- string
405 function o = getBoundaryOperator(obj, op, boundary)
406 assertIsMember(boundary, {'w', 'e', 's', 'n'})
407 assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
408
409 o = obj.([op, '_', boundary]);
410
411 end
412
413 % Returns the boundary operator op for the boundary specified by the string boundary.
414 % op -- string
415 function o = getBoundaryOperatorForScalarField(obj, op, boundary)
416 assertIsMember(boundary, {'w', 'e', 's', 'n'})
417 assertIsMember(op, {'e'})
418
419 switch op
420
421 case 'e'
422 o = obj.(['e_scalar', '_', boundary]);
423 end
424
425 end
426
427 % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary.
428 % Formula: tau_i = T_ij u_j
429 % op -- string
430 function T = getBoundaryTractionOperator(obj, boundary)
431 assertIsMember(boundary, {'w', 'e', 's', 'n'})
432
433 T = obj.(['T', '_', boundary]);
434 end
435
436 % Returns square boundary quadrature matrix, of dimension
437 % corresponding to the number of boundary unknowns
438 %
439 % boundary -- string
440 function H = getBoundaryQuadrature(obj, boundary)
441 assertIsMember(boundary, {'w', 'e', 's', 'n'})
442
443 H = obj.getBoundaryQuadratureForScalarField(boundary);
444 I_dim = speye(obj.dim, obj.dim);
445 H = kron(H, I_dim);
446 end
447
448 % Returns square boundary quadrature matrix, of dimension
449 % corresponding to the number of boundary grid points
450 %
451 % boundary -- string
452 function H_b = getBoundaryQuadratureForScalarField(obj, boundary)
453 assertIsMember(boundary, {'w', 'e', 's', 'n'})
454
455 H_b = obj.(['H_', boundary]);
456 end
457
458 function N = size(obj)
459 N = obj.dim*prod(obj.m);
460 end
461 end
462 end