Mercurial > repos > public > sbplib

comparison +scheme/elasticShearVariable.m @ 663:b45ec2b28cc2 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
First implementation of elastic shear operator with free boundary BC.
author Martin Almquist <malmquist@stanford.edu>
date Thu, 14 Dec 2017 13:54:20 -0800
parents
children 8e6dfd22fc59
comparison
equal deleted inserted replaced
662:b189bc409cdb 663:b45ec2b28cc2
1 classdef elasticShearVariable < scheme.Scheme
2 properties
3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing
5
6 grid
7
8 order % Order accuracy for the approximation
9
10 a % Variable coefficient lambda of the operator
11 rho % Density
12
13 D % Total operator
14 D1 % First derivatives
15 D2 % Second derivatives
16 H, Hi % Inner products
17 e_l, e_r
18 d_l, d_r % Normal derivatives at the boundary
19 H_boundary % Boundary inner products
20 end
21
22 methods
23 % Implements the shear part of the elastic wave equation, i.e.
24 % rho u_{i,tt} = d_i a d_j u_j + d_j a d_j u_i
25 % where a = lambda.
26
27 function obj = elasticShearVariable(g ,order, a_fun, rho_fun, opSet)
28 default_arg('opSet',@sbp.D2Variable);
29 default_arg('a_fun', @(x,y) 0*x+1);
30 default_arg('rho_fun', @(x,y) 0*x+1);
31 dim = 2;
32
33 assert(isa(g, 'grid.Cartesian'))
34
35 a = grid.evalOn(g, a_fun);
36 rho = grid.evalOn(g, rho_fun);
37 m = g.size();
38 m_tot = g.N();
39
40 h = g.scaling();
41
42 % 1D operators
43 ops = cell(dim,1);
44 for i = 1:dim
45 ops{i} = opSet(m(i), {0, 1}, order);
46 end
47
48 I = cell(dim,1);
49 D1 = cell(dim,1);
50 D2 = cell(dim,1);
51 H = cell(dim,1);
52 Hi = cell(dim,1);
53 e_l = cell(dim,1);
54 e_r = cell(dim,1);
55 d1_l = cell(dim,1);
56 d1_r = cell(dim,1);
57
58 for i = 1:dim
59 I{i} = speye{m(i));
60 D1{i} = ops{i}.D1;
61 D2{i} = ops{i}.D2;
62 H{i} = ops{i}.H;
63 Hi{i} = ops{i}.HI;
64 e_l{i} = ops{i}.e_l;
65 e_r{i} = ops{i}.e_r;
66 d1_l{i} = ops{i}.d1_l;
67 d1_r{i} = ops{i}.d1_r;
68 end
69
70 %====== Assemble full operators ========
71 A = spdiag(a);
72 RHO = spdiag(rho);
73
74 obj.D1 = cell(dim,1);
75 obj.D2 = cell(dim,1);
76 obj.e_l = cell(dim,1);
77 obj.e_r = cell(dim,1);
78 obj.d1_l = cell(dim,1);
79 obj.d1_r = cell(dim,1);
80
81 % D1
82 obj.D1{1} = kron{D1{1},I(2)};
83 obj.D1{2} = kron{I{1},D1(2)};
84
85 % Boundary operators
86 obj.e_l{1} = kron{e_l{1},I(2)};
87 obj.e_l{2} = kron{I{1},e_l(2)};
88 obj.e_r{1} = kron{e_r{1},I(2)};
89 obj.e_r{2} = kron{I{1},e_r(2)};
90
91 obj.d1_l{1} = kron{d1_l{1},I(2)};
92 obj.d1_l{2} = kron{I{1},d1_l(2)};
93 obj.d1_r{1} = kron{d1_r{1},I(2)};
94 obj.d1_r{2} = kron{I{1},d1_r(2)};
95
96 % D2
97 for i = 1:dim
98 obj.D2{i} = sparse(m_tot);
99 end
100 ind = grid.funcToMatrix(g, 1:m_tot);
101
102 for i = 1:m(2)
103 D = D2{1}(a(ind(:,i)));
104 p = ind(:,i);
105 obj.D2{1}(p,p) = D;
106 end
107
108 for i = 1:m(1)
109 D = D2{2}(a(ind(i,:)));
110 p = ind(i,:);
111 obj.D2{2}(p,p) = D;
112 end
113
114 % Quadratures
115 obj.H = kron(H{1},H{2});
116 obj.H_boundary = cell(dim,1);
117 obj.H_boundary{1} = H{2};
118 obj.H_boundary{2} = H{1};
119
120 % Boundary coefficient matrices and quadratures
121 obj.A_boundary_l = cell(dim,1);
122 obj.A_boundary_r = cell(dim,1);
123 for i = 1:dim
124 obj.A_boundary_l{i} = e_l{i}'*A*e_l{i};
125 obj.A_boundary_r{i} = e_r{i}'*A*e_r{i};
126 end
127
128 % E{i}^T picks out component i.
129 E = cell(dim,1);
130 I = speye{mtot,mtot};
131 for i = 1:dim
132 e = sparse(dim,1);
133 e(i) = 1;
134 E{i} = kron(e,I)
135 end
136 obj.E = E;
137
138 % Differentiation matrix D (without SAT)
139 D = 0;
140 d = @kroneckerDelta; % Kronecker delta
141 db = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta
142 for i = 1:dim
143 for j = 1:dim
144 D = D + E{i}*inv(rho)*( d(i,j)*D2{i}*E{j}' +...
145 db(i,j)*D1{j}*A*D1{i}*E{j}' + ...
146 D2{j}*E{i}' ...
147 );
148 end
149 end
150 obj.D = D;
151 %=========================================%
152
153
154
155 % Misc.
156 obj.m = m;
157 obj.h = h;
158 obj.order = order;
159 obj.grid = g;
160
161 obj.a = a;
162 obj.b = b;
163
164 % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
165 % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
166 end
167
168
169 % Closure functions return the operators applied to the own domain to close the boundary
170 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
171 % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
172 % type is a string specifying the type of boundary condition if there are several.
173 % data is a function returning the data that should be applied at the boundary.
174 % neighbour_scheme is an instance of Scheme that should be interfaced to.
175 % neighbour_boundary is a string specifying which boundary to interface to.
176 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
177 default_arg('type','free');
178 default_arg('parameter', []);
179
180 delta = @kroneckerDelta; % Kronecker delta
181 delta_b = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta
182
183 % j is the coordinate direction of the boundary
184 % nj: outward unit normal component.
185 % nj = -1 for west, south, bottom boundaries
186 % nj = 1 for east, north, top boundaries
187 [j, nj] = obj.get_boundary_number(boundary);
188 switch nj
189 case 1
190 e = obj.e_r;
191 d = obj.d_r;
192 case -1
193 e = obj.e_l;
194 d = obj.d_l;
195 end
196
197 E = obj.E;
198 Hi = obj.Hi;
199 H_gamma = obj.H_boundary{j};
200 A = obj.A;
201
202 switch type
203 % Dirichlet boundary condition
204 case {'D','d','dirichlet'}
205 error('Dirichlet not implemented')
206 tuning = 1.2;
207 % tuning = 20.2;
208
209 b1 = gamm*obj.lambda./obj.a11.^2;
210 b2 = gamm*obj.lambda./obj.a22.^2;
211
212 tau1 = tuning * spdiag(-1./b1 - 1./b2);
213 tau2 = 1;
214
215 tau = (tau1*e + tau2*d)*H_b;
216
217 closure = obj.a*obj.Hi*tau*e';
218 penalty = -obj.a*obj.Hi*tau;
219
220
221 % Free boundary condition
222 case {'F','f','Free','free'}
223 closure = 0;
224 penalty = 0;
225 % Loop over components
226 for i = 1:3
227 closure = closure + E{i}*(-sign)*Hi*e{j}*H_gamma*(...
228 e{j}'*A*e{j}*d{j}'*E{i}' + ...
229 delta(i,j)*e{j}'*A*e{i}*d{i}*E{j}' + ...
230 delta_b(i,j)*A*D1{i}*E{j}' ...
231 );
232 penalty = penalty - E{i}*(-sign)*Hi*e{j}*H_gamma;
233 end
234
235 % Unknown boundary condition
236 otherwise
237 error('No such boundary condition: type = %s',type);
238 end
239 end
240
241 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
242 % u denotes the solution in the own domain
243 % v denotes the solution in the neighbour domain
244 tuning = 1.2;
245 % tuning = 20.2;
246 error('Interface not implemented');
247 end
248
249 % Ruturns the coordinate number and outward normal component for the boundary specified by the string boundary.
250 function [j, nj] = get_boundary_number(obj, boundary)
251
252 switch boundary
253 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
254 j = 1;
255 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
256 j = 2;
257 otherwise
258 error('No such boundary: boundary = %s',boundary);
259 end
260
261 switch boundary
262 case {'w','W','west','West','s','S','south','South'}
263 nj = -1;
264 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
265 nj = 1;
266 end
267 end
268
269 function N = size(obj)
270 N = prod(obj.m);
271 end
272 end
273 end