Mercurial > repos > public > sbplib

comparison +scheme/elasticShearVariable.m @ 669:17e62551cdc2 feature/poroelastic

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add Dirichlet condition. Symmetric, semi-definite matrix and MMS convergence working with variable coeff!
author Martin Almquist <malmquist@stanford.edu>
date Fri, 22 Dec 2017 14:16:17 +0100
parents ed853945ee99
children 9e1d2351f539
comparison
equal deleted inserted replaced
668:43ea848b6aa1 669:17e62551cdc2
13 13
14 D % Total operator 14 D % Total operator
15 D1 % First derivatives 15 D1 % First derivatives
16 D2 % Second derivatives 16 D2 % Second derivatives
17 H, Hi % Inner products 17 H, Hi % Inner products
18 phi % Borrowing constant for (d1 - e^T*D1) from R
19 gamma % Borrowing constant for d1 from M
20 H11 % First element of H
18 e_l, e_r 21 e_l, e_r
19 d1_l, d1_r % Normal derivatives at the boundary 22 d1_l, d1_r % Normal derivatives at the boundary
20 E % E{i}^T picks out component i 23 E % E{i}^T picks out component i
21 24
22 H_boundary % Boundary inner products 25 H_boundary % Boundary inner products
49 % 1D operators 52 % 1D operators
50 ops = cell(dim,1); 53 ops = cell(dim,1);
51 for i = 1:dim 54 for i = 1:dim
52 ops{i} = opSet(m(i), {0, L(i)}, order); 55 ops{i} = opSet(m(i), {0, L(i)}, order);
53 end 56 end
57
58 % Borrowing constants
59 beta = ops{1}.borrowing.R.delta_D;
60 obj.H11 = ops{1}.borrowing.H11;
61 obj.phi = beta/obj.H11;
62 obj.gamma = ops{1}.borrowing.M.d1;
54 63
55 I = cell(dim,1); 64 I = cell(dim,1);
56 D1 = cell(dim,1); 65 D1 = cell(dim,1);
57 D2 = cell(dim,1); 66 D2 = cell(dim,1);
58 H = cell(dim,1); 67 H = cell(dim,1);
208 H_gamma = obj.H_boundary{j}; 217 H_gamma = obj.H_boundary{j};
209 A = obj.A; 218 A = obj.A;
210 RHO = obj.RHO; 219 RHO = obj.RHO;
211 D1 = obj.D1; 220 D1 = obj.D1;
212 221
222 phi = obj.phi;
223 gamma = obj.gamma;
224 H11 = obj.H11;
225 h = obj.h;
226
213 switch type 227 switch type
214 % Dirichlet boundary condition 228 % Dirichlet boundary condition
215 case {'D','d','dirichlet'} 229 case {'D','d','dirichlet'}
216 error('Dirichlet not implemented')
217 tuning = 1.2; 230 tuning = 1.2;
218 % tuning = 20.2; 231 phi = obj.phi;
219 232
220 b1 = gamm*obj.lambda./obj.a11.^2; 233 closures = cell(obj.dim,1);
221 b2 = gamm*obj.lambda./obj.a22.^2; 234 penalties = cell(obj.dim,obj.dim);
222 235 % Loop over components
223 tau1 = tuning * spdiag(-1./b1 - 1./b2); 236 for i = 1:obj.dim
224 tau2 = 1; 237 H_gamma_i = obj.H_boundary{i};
225 238 sigma_ij = tuning*delta(i,j)*2/(gamma*h(j)) +...
226 tau = (tau1*e + tau2*d)*H_b; 239 tuning*delta_b(i,j)*(2/(H11*h(j)) + 2/(H11*h(j)*phi));
227 240
228 closure = obj.a*obj.Hi*tau*e'; 241 ci = E{i}*inv(RHO)*nj*Hi*...
229 penalty = -obj.a*obj.Hi*tau; 242 ( (e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{i}' + ...
230 243 delta(i,j)*(e{j}*H_gamma*e{j}'*A*e{j}*d{j}')'*E{j}' ...
244 ) ...
245 - sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma*e{j}'*E{i}';
246
247 cj = E{j}*inv(RHO)*nj*Hi*...
248 ( delta_b(i,j)*(e{j}*H_gamma*e{j}'*A*D1{i})'*E{i}' ...
249 );
250
251 if isempty(closures{i})
252 closures{i} = ci;
253 else
254 closures{i} = closures{i} + ci;
255 end
256
257 if isempty(closures{j})
258 closures{j} = cj;
259 else
260 closures{j} = closures{j} + cj;
261 end
262
263 pii = - E{i}*inv(RHO)*nj*Hi*...
264 ( (H_gamma*e{j}'*A*e{j}*d{j}')' + ...
265 delta(i,j)*(H_gamma*e{j}'*A*e{j}*d{j}')' ...
266 ) ...
267 + sigma_ij*E{i}*inv(RHO)*Hi*A*e{j}*H_gamma;
268
269 pji = - E{j}*inv(RHO)*nj*Hi*...
270 ( delta_b(i,j)*(H_gamma*e{j}'*A*D1{i})' );
271
272 % Dummies
273 pij = - 0*E{i}*e{j};
274 pjj = - 0*E{j}*e{j};
275
276 if isempty(penalties{i,i})
277 penalties{i,i} = pii;
278 else
279 penalties{i,i} = penalties{i,i} + pii;
280 end
281
282 if isempty(penalties{j,i})
283 penalties{j,i} = pji;
284 else
285 penalties{j,i} = penalties{j,i} + pji;
286 end
287
288 if isempty(penalties{i,j})
289 penalties{i,j} = pij;
290 else
291 penalties{i,j} = penalties{i,j} + pij;
292 end
293
294 if isempty(penalties{j,j})
295 penalties{j,j} = pjj;
296 else
297 penalties{j,j} = penalties{j,j} + pjj;
298 end
299 end
300 [rows, cols] = size(closures{1});
301 closure = sparse(rows, cols);
302 for i = 1:obj.dim
303 closure = closure + closures{i};
304 end
305 penalty = penalties;
231 306
232 % Free boundary condition 307 % Free boundary condition
233 case {'F','f','Free','free'} 308 case {'F','f','Free','free'}
234 closures = cell(obj.dim,1); 309 closures = cell(obj.dim,1);
235 penalties = cell(obj.dim,1); 310 penalties = cell(obj.dim,1);