Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dCurvilinearAnisotropic.m @ 1205:3258dca12af8 feature/poroelastic

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