Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1086:90c23fd08f59 feature/laplace_curvilinear_test

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make LaplaceCurvilinearNew the default version and remove the others
author Martin Almquist <malmquist@stanford.edu>
date Fri, 29 Mar 2019 14:41:36 -0700
parents 1341c6cea9c1
children 867307f4d80f
comparison
equal deleted inserted replaced
1085:49c0b8c7330a 1086:90c23fd08f59
1 classdef LaplaceCurvilinear < scheme.Scheme 1 classdef LaplaceCurvilinear < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 dim % Number of spatial dimensions
5 6
6 grid 7 grid
7 8
8 order % Order accuracy for the approximation 9 order % Order of accuracy for the approximation
9 10
10 a,b % Parameters of the operator 11 a,b % Parameters of the operator
11 12
12 13
13 % Inner products and operators for physical coordinates 14 % Inner products and operators for physical coordinates
20 M % Gradient inner product 21 M % Gradient inner product
21 22
22 % Metric coefficients 23 % Metric coefficients
23 J, Ji 24 J, Ji
24 a11, a12, a22 25 a11, a12, a22
26 K
25 x_u 27 x_u
26 x_v 28 x_v
27 y_u 29 y_u
28 y_v 30 y_v
31 s_w, s_e, s_s, s_n % Boundary integral scale factors
29 32
30 % Inner product and operators for logical coordinates 33 % Inner product and operators for logical coordinates
31 H_u, H_v % Norms in the x and y directions 34 H_u, H_v % Norms in the x and y directions
32 Hi_u, Hi_v 35 Hi_u, Hi_v
33 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. 36 Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir.
34 Hiu, Hiv 37 Hiu, Hiv
35 du_w, dv_w 38 du_w, dv_w
36 du_e, dv_e 39 du_e, dv_e
37 du_s, dv_s 40 du_s, dv_s
38 du_n, dv_n 41 du_n, dv_n
42
43 % Borrowing constants
44 theta_M_u, theta_M_v
45 theta_R_u, theta_R_v
46 theta_H_u, theta_H_v
47
48 % Temporary
49 lambda
39 gamm_u, gamm_v 50 gamm_u, gamm_v
40 lambda
41 51
42 end 52 end
43 53
44 methods 54 methods
45 % Implements a*div(b*grad(u)) as a SBP scheme 55 % Implements a*div(b*grad(u)) as a SBP scheme
46 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?) 56 % TODO: Implement proper H, it should be the real physical quadrature, the logic quadrature may be but in a separate variable (H_logic?)
47 57
48 function obj = LaplaceCurvilinear(g ,order, a, b, opSet) 58 function obj = LaplaceCurvilinear(g, order, a, b, opSet)
49 default_arg('opSet',@sbp.D2Variable); 59 default_arg('opSet',@sbp.D2Variable);
50 default_arg('a', 1); 60 default_arg('a', 1);
51 default_arg('b', 1); 61 default_arg('b', @(x,y) 0*x + 1);
52
53 if b ~=1
54 error('Not implemented yet')
55 end
56 62
57 % assert(isa(g, 'grid.Curvilinear')) 63 % assert(isa(g, 'grid.Curvilinear'))
58 if isa(a, 'function_handle') 64 if isa(a, 'function_handle')
59 a = grid.evalOn(g, a); 65 a = grid.evalOn(g, a);
60 a = spdiag(a); 66 a = spdiag(a);
61 end 67 end
62 68
69 if isa(b, 'function_handle')
70 b = grid.evalOn(g, b);
71 b = spdiag(b);
72 end
73
74 % If b is scalar
75 if length(b) == 1
76 b = b*speye(g.N(), g.N());
77 end
78
79 dim = 2;
63 m = g.size(); 80 m = g.size();
64 m_u = m(1); 81 m_u = m(1);
65 m_v = m(2); 82 m_v = m(2);
66 m_tot = g.N(); 83 m_tot = g.N();
67 84
132 a11 = 1./J .* (x_v.^2 + y_v.^2); 149 a11 = 1./J .* (x_v.^2 + y_v.^2);
133 a12 = -1./J .* (x_u.*x_v + y_u.*y_v); 150 a12 = -1./J .* (x_u.*x_v + y_u.*y_v);
134 a22 = 1./J .* (x_u.^2 + y_u.^2); 151 a22 = 1./J .* (x_u.^2 + y_u.^2);
135 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); 152 lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2));
136 153
154 K = cell(dim, dim);
155 K{1,1} = spdiag(y_v./J);
156 K{1,2} = spdiag(-y_u./J);
157 K{2,1} = spdiag(-x_v./J);
158 K{2,2} = spdiag(x_u./J);
159 obj.K = K;
160
137 obj.x_u = x_u; 161 obj.x_u = x_u;
138 obj.x_v = x_v; 162 obj.x_v = x_v;
139 obj.y_u = y_u; 163 obj.y_u = y_u;
140 obj.y_v = y_v; 164 obj.y_v = y_v;
141 165
142
143 % Assemble full operators 166 % Assemble full operators
144 L_12 = spdiag(a12); 167 L_12 = spdiag(a12);
145 Duv = Du*L_12*Dv; 168 Duv = Du*b*L_12*Dv;
146 Dvu = Dv*L_12*Du; 169 Dvu = Dv*b*L_12*Du;
147 170
148 Duu = sparse(m_tot); 171 Duu = sparse(m_tot);
149 Dvv = sparse(m_tot); 172 Dvv = sparse(m_tot);
150 ind = grid.funcToMatrix(g, 1:m_tot); 173 ind = grid.funcToMatrix(g, 1:m_tot);
151 174
152 for i = 1:m_v 175 for i = 1:m_v
153 D = D2_u(a11(ind(:,i))); 176 b_a11 = b*a11;
177 D = D2_u(b_a11(ind(:,i)));
154 p = ind(:,i); 178 p = ind(:,i);
155 Duu(p,p) = D; 179 Duu(p,p) = D;
156 end 180 end
157 181
158 for i = 1:m_u 182 for i = 1:m_u
159 D = D2_v(a22(ind(i,:))); 183 b_a22 = b*a22;
184 D = D2_v(b_a22(ind(i,:)));
160 p = ind(i,:); 185 p = ind(i,:);
161 Dvv(p,p) = D; 186 Dvv(p,p) = D;
162 end 187 end
163 188
164 189
212 % Misc. 237 % Misc.
213 obj.m = m; 238 obj.m = m;
214 obj.h = [h_u h_v]; 239 obj.h = [h_u h_v];
215 obj.order = order; 240 obj.order = order;
216 obj.grid = g; 241 obj.grid = g;
242 obj.dim = dim;
217 243
218 obj.a = a; 244 obj.a = a;
219 obj.b = b; 245 obj.b = b;
220 obj.a11 = a11; 246 obj.a11 = a11;
221 obj.a12 = a12; 247 obj.a12 = a12;
222 obj.a22 = a22; 248 obj.a22 = a22;
249 obj.s_w = spdiag(s_w);
250 obj.s_e = spdiag(s_e);
251 obj.s_s = spdiag(s_s);
252 obj.s_n = spdiag(s_n);
253
254 obj.theta_M_u = h_u*ops_u.borrowing.M.d1;
255 obj.theta_M_v = h_v*ops_v.borrowing.M.d1;
256
257 obj.theta_R_u = h_u*ops_u.borrowing.R.delta_D;
258 obj.theta_R_v = h_v*ops_v.borrowing.R.delta_D;
259
260 obj.theta_H_u = h_u*ops_u.borrowing.H11;
261 obj.theta_H_v = h_v*ops_v.borrowing.H11;
262
263 % Temporary
223 obj.lambda = lambda; 264 obj.lambda = lambda;
224
225 obj.gamm_u = h_u*ops_u.borrowing.M.d1; 265 obj.gamm_u = h_u*ops_u.borrowing.M.d1;
226 obj.gamm_v = h_v*ops_v.borrowing.M.d1; 266 obj.gamm_v = h_v*ops_v.borrowing.M.d1;
227 end 267 end
228 268
229 269
236 % neighbour_boundary is a string specifying which boundary to interface to. 276 % neighbour_boundary is a string specifying which boundary to interface to.
237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 277 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
238 default_arg('type','neumann'); 278 default_arg('type','neumann');
239 default_arg('parameter', []); 279 default_arg('parameter', []);
240 280
241 [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); 281 e = obj.getBoundaryOperator('e', boundary);
282 d = obj.getBoundaryOperator('d', boundary);
242 H_b = obj.getBoundaryQuadrature(boundary); 283 H_b = obj.getBoundaryQuadrature(boundary);
243 gamm = obj.getBoundaryBorrowing(boundary); 284 s_b = obj.getBoundaryScaling(boundary);
285 [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary);
286 m = obj.getBoundaryNumber(boundary);
287
288 K = obj.K;
289 J = obj.J;
290 Hi = obj.Hi;
291 a = obj.a;
292 b_b = e'*obj.b*e;
293
244 switch type 294 switch type
245 % Dirichlet boundary condition 295 % Dirichlet boundary condition
246 case {'D','d','dirichlet'} 296 case {'D','d','dirichlet'}
247 tuning = 1.0; 297 tuning = 1.0;
248 298
249 b1 = gamm*obj.lambda./obj.a11.^2; 299 sigma = 0*b_b;
250 b2 = gamm*obj.lambda./obj.a22.^2; 300 for i = 1:obj.dim
251 301 sigma = sigma + e'*J*K{i,m}*K{i,m}*e;
252 tau1 = tuning * spdiag(-1./b1 - 1./b2); 302 end
253 tau2 = 1; 303 sigma = sigma/s_b;
254 304 tau = tuning*(1/th_R + obj.dim/th_H)*sigma;
255 tau = (tau1*e + tau2*d)*H_b; 305
256 306 closure = a*Hi*d*b_b*H_b*e' ...
257 closure = obj.a*obj.Hi*tau*e'; 307 -a*Hi*e*tau*b_b*H_b*e';
258 penalty = -obj.a*obj.Hi*tau; 308
259 309 penalty = -a*Hi*d*b_b*H_b ...
260 310 +a*Hi*e*tau*b_b*H_b;
261 % Neumann boundary condition 311
312
313 % Neumann boundary condition. Note that the penalty is for du/dn and not b*du/dn.
262 case {'N','n','neumann'} 314 case {'N','n','neumann'}
263 tau1 = -1; 315 tau1 = -1;
264 tau2 = 0; 316 tau2 = 0;
265 tau = (tau1*e + tau2*d)*H_b; 317 tau = (tau1*e + tau2*d)*H_b;
266 318
267 closure = obj.a*obj.Hi*tau*d'; 319 closure = a*Hi*tau*b_b*d';
268 penalty = -obj.a*obj.Hi*tau; 320 penalty = -a*Hi*tau*b_b;
269 321
270 322
271 % Unknown, boundary condition 323 % Unknown, boundary condition
272 otherwise 324 otherwise
273 error('No such boundary condition: type = %s',type); 325 error('No such boundary condition: type = %s',type);
278 % Fields: 330 % Fields:
279 % -- tuning: penalty strength, defaults to 1.2 331 % -- tuning: penalty strength, defaults to 1.2
280 % -- interpolation: type of interpolation, default 'none' 332 % -- interpolation: type of interpolation, default 'none'
281 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) 333 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type)
282 334
283 defaultType.tuning = 1.2; 335 % error('Not implemented')
336
337 defaultType.tuning = 1.0;
284 defaultType.interpolation = 'none'; 338 defaultType.interpolation = 'none';
285 default_struct('type', defaultType); 339 default_struct('type', defaultType);
286 340
287 switch type.interpolation 341 switch type.interpolation
288 case {'none', ''} 342 case {'none', ''}
295 end 349 end
296 350
297 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 351 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
298 tuning = type.tuning; 352 tuning = type.tuning;
299 353
354 dim = obj.dim;
300 % u denotes the solution in the own domain 355 % u denotes the solution in the own domain
301 % v denotes the solution in the neighbour domain 356 % v denotes the solution in the neighbour domain
302 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); 357 u = obj;
358 v = neighbour_scheme;
359
360 % Boundary operators, u
361 e_u = u.getBoundaryOperator('e', boundary);
362 d_u = u.getBoundaryOperator('d', boundary);
363 gamm_u = u.getBoundaryBorrowing(boundary);
364 s_b_u = u.getBoundaryScaling(boundary);
365 [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary);
366 m_u = u.getBoundaryNumber(boundary);
367
368 % Coefficients, u
369 K_u = u.K;
370 J_u = u.J;
371 b_b_u = e_u'*u.b*e_u;
372
373 % Boundary operators, v
374 e_v = v.getBoundaryOperator('e', neighbour_boundary);
375 d_v = v.getBoundaryOperator('d', neighbour_boundary);
376 gamm_v = v.getBoundaryBorrowing(neighbour_boundary);
377 s_b_v = v.getBoundaryScaling(neighbour_boundary);
378 [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary);
379 m_v = v.getBoundaryNumber(neighbour_boundary);
380
381 % Coefficients, v
382 K_v = v.K;
383 J_v = v.J;
384 b_b_v = e_v'*v.b*e_v;
385
386 %--- Penalty strength tau -------------
387 sigma_u = 0*b_b_u;
388 sigma_v = 0*b_b_v;
389 for i = 1:obj.dim
390 sigma_u = sigma_u + e_u'*J_u*K_u{i,m_u}*K_u{i,m_u}*e_u;
391 sigma_v = sigma_v + e_v'*J_v*K_v{i,m_v}*K_v{i,m_v}*e_v;
392 end
393 sigma_u = sigma_u/s_b_u;
394 sigma_v = sigma_v/s_b_v;
395
396 tau_R_u = 1/th_R_u*sigma_u;
397 tau_R_v = 1/th_R_v*sigma_v;
398
399 tau_H_u = dim*1/th_H_u*sigma_u;
400 tau_H_v = dim*1/th_H_v*sigma_v;
401
402 tau = 1/4*tuning*(b_b_u*(tau_R_u + tau_H_u) + b_b_v*(tau_R_v + tau_H_v));
403 %--------------------------------------
404
405 % Operators/coefficients that are only required from this side
406 Hi = u.Hi;
407 H_b = u.getBoundaryQuadrature(boundary);
408 a = u.a;
409
410 closure = 1/2*a*Hi*d_u*b_b_u*H_b*e_u' ...
411 -1/2*a*Hi*e_u*H_b*b_b_u*d_u' ...
412 -a*Hi*e_u*tau*H_b*e_u';
413
414 penalty = -1/2*a*Hi*d_u*b_b_u*H_b*e_v' ...
415 -1/2*a*Hi*e_u*H_b*b_b_v*d_v' ...
416 +a*Hi*e_u*tau*H_b*e_v';
417 end
418
419 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
420
421 % TODO: Make this work for curvilinear grids
422 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
423 warning('LaplaceCurvilinear: Non-conforming interface uses Virtas penalty strength');
424 warning('LaplaceCurvilinear: Non-conforming interface assumes that b is constant');
425
426 % User can request special interpolation operators by specifying type.interpOpSet
427 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
428 interpOpSet = type.interpOpSet;
429 tuning = type.tuning;
430
431
432 % u denotes the solution in the own domain
433 % v denotes the solution in the neighbour domain
434 e_u = obj.getBoundaryOperator('e', boundary);
435 d_u = obj.getBoundaryOperator('d', boundary);
303 H_b_u = obj.getBoundaryQuadrature(boundary); 436 H_b_u = obj.getBoundaryQuadrature(boundary);
304 I_u = obj.getBoundaryIndices(boundary); 437 I_u = obj.getBoundaryIndices(boundary);
305 gamm_u = obj.getBoundaryBorrowing(boundary); 438 gamm_u = obj.getBoundaryBorrowing(boundary);
306 439
307 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); 440 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
441 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
308 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); 442 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
309 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); 443 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
310 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); 444 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
445
446
447 % Find the number of grid points along the interface
448 m_u = size(e_u, 2);
449 m_v = size(e_v, 2);
450
451 Hi = obj.Hi;
452 a = obj.a;
311 453
312 u = obj; 454 u = obj;
313 v = neighbour_scheme; 455 v = neighbour_scheme;
314 456
315 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 457 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
316 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; 458 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
317 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; 459 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
318 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; 460 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
319 461
320 tau1 = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v);
321 tau1 = tuning * spdiag(tau1);
322 tau2 = 1/2;
323
324 sig1 = -1/2;
325 sig2 = 0;
326
327 tau = (e_u*tau1 + tau2*d_u)*H_b_u;
328 sig = (sig1*e_u + sig2*d_u)*H_b_u;
329
330 closure = obj.a*obj.Hi*( tau*e_u' + sig*d_u');
331 penalty = obj.a*obj.Hi*(-tau*e_v' + sig*d_v');
332 end
333
334 function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
335
336 % TODO: Make this work for curvilinear grids
337 warning('LaplaceCurvilinear: Non-conforming grid interpolation only works for Cartesian grids.');
338
339 % User can request special interpolation operators by specifying type.interpOpSet
340 default_field(type, 'interpOpSet', @sbp.InterpOpsOP);
341 interpOpSet = type.interpOpSet;
342 tuning = type.tuning;
343
344
345 % u denotes the solution in the own domain
346 % v denotes the solution in the neighbour domain
347 [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary);
348 H_b_u = obj.getBoundaryQuadrature(boundary);
349 I_u = obj.getBoundaryIndices(boundary);
350 gamm_u = obj.getBoundaryBorrowing(boundary);
351
352 [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary);
353 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
354 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
355 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
356
357
358 % Find the number of grid points along the interface
359 m_u = size(e_u, 2);
360 m_v = size(e_v, 2);
361
362 Hi = obj.Hi;
363 a = obj.a;
364
365 u = obj;
366 v = neighbour_scheme;
367
368 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
369 b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2;
370 b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2;
371 b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2;
372
373 tau_u = -1./(4*b1_u) -1./(4*b2_u); 462 tau_u = -1./(4*b1_u) -1./(4*b2_u);
374 tau_v = -1./(4*b1_v) -1./(4*b2_v); 463 tau_v = -1./(4*b1_v) -1./(4*b2_v);
375 464
376 tau_u = tuning * spdiag(tau_u); 465 tau_u = tuning * spdiag(tau_u);
377 tau_v = tuning * spdiag(tau_v); 466 tau_v = tuning * spdiag(tau_v);
393 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; 482 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
394 483
395 end 484 end
396 485
397 % Returns the boundary operator op for the boundary specified by the string boundary. 486 % Returns the boundary operator op for the boundary specified by the string boundary.
398 % op -- string or a cell array of strings 487 % op -- string
399 % boundary -- string 488 % boundary -- string
400 function varargout = getBoundaryOperator(obj, op, boundary) 489 function o = getBoundaryOperator(obj, op, boundary)
401 490 assertIsMember(op, {'e', 'd'})
402 if ~iscell(op) 491 assertIsMember(boundary, {'w', 'e', 's', 'n'})
403 op = {op}; 492
404 end 493 o = obj.([op, '_', boundary]);
405
406 for i = 1:numel(op)
407 switch op{i}
408 case 'e'
409 switch boundary
410 case 'w'
411 e = obj.e_w;
412 case 'e'
413 e = obj.e_e;
414 case 's'
415 e = obj.e_s;
416 case 'n'
417 e = obj.e_n;
418 otherwise
419 error('No such boundary: boundary = %s',boundary);
420 end
421 varargout{i} = e;
422
423 case 'd'
424 switch boundary
425 case 'w'
426 d = obj.d_w;
427 case 'e'
428 d = obj.d_e;
429 case 's'
430 d = obj.d_s;
431 case 'n'
432 d = obj.d_n;
433 otherwise
434 error('No such boundary: boundary = %s',boundary);
435 end
436 varargout{i} = d;
437 end
438 end
439 end 494 end
440 495
441 % Returns square boundary quadrature matrix, of dimension 496 % Returns square boundary quadrature matrix, of dimension
442 % corresponding to the number of boundary points 497 % corresponding to the number of boundary points
443 % 498 %
444 % boundary -- string 499 % boundary -- string
445 function H_b = getBoundaryQuadrature(obj, boundary) 500 function H_b = getBoundaryQuadrature(obj, boundary)
501 assertIsMember(boundary, {'w', 'e', 's', 'n'})
502
503 H_b = obj.(['H_', boundary]);
504 end
505
506 % Returns square boundary quadrature scaling matrix, of dimension
507 % corresponding to the number of boundary points
508 %
509 % boundary -- string
510 function s_b = getBoundaryScaling(obj, boundary)
511 assertIsMember(boundary, {'w', 'e', 's', 'n'})
512
513 s_b = obj.(['s_', boundary]);
514 end
515
516 % Returns the coordinate number corresponding to the boundary
517 %
518 % boundary -- string
519 function m = getBoundaryNumber(obj, boundary)
520 assertIsMember(boundary, {'w', 'e', 's', 'n'})
446 521
447 switch boundary 522 switch boundary
448 case 'w' 523 case {'w', 'e'}
449 H_b = obj.H_w; 524 m = 1;
450 case 'e' 525 case {'s', 'n'}
451 H_b = obj.H_e; 526 m = 2;
452 case 's'
453 H_b = obj.H_s;
454 case 'n'
455 H_b = obj.H_n;
456 otherwise
457 error('No such boundary: boundary = %s',boundary);
458 end 527 end
459 end 528 end
460 529
461 % Returns the indices of the boundary points in the grid matrix 530 % Returns the indices of the boundary points in the grid matrix
462 % boundary -- string 531 % boundary -- string
476 end 545 end
477 end 546 end
478 547
479 % Returns borrowing constant gamma 548 % Returns borrowing constant gamma
480 % boundary -- string 549 % boundary -- string
481 function gamm = getBoundaryBorrowing(obj, boundary) 550 function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary)
482 switch boundary 551 switch boundary
483 case {'w','e'} 552 case {'w','e'}
484 gamm = obj.gamm_u; 553 theta_H = obj.theta_H_u;
554 theta_M = obj.theta_M_u;
555 theta_R = obj.theta_R_u;
485 case {'s','n'} 556 case {'s','n'}
486 gamm = obj.gamm_v; 557 theta_H = obj.theta_H_v;
558 theta_M = obj.theta_M_v;
559 theta_R = obj.theta_R_v;
487 otherwise 560 otherwise
488 error('No such boundary: boundary = %s',boundary); 561 error('No such boundary: boundary = %s',boundary);
489 end 562 end
490 end 563 end
491 564