Mercurial > repos > public > sbplib

comparison +scheme/Heat2dCurvilinear.m @ 997:78db023a7fe3 feature/getBoundaryOp

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add getBoundaryOperator method to all 2d schemes, except obsolete Wave2dCurve and ElastiCurve, which needs a big makeover.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 12 Jan 2019 11:57:50 -0800
parents 08f3ffe63f48
children 8d73fcdb07a5
comparison
equal deleted inserted replaced
982:a0b3161e44f3 997:78db023a7fe3
1 classdef Heat2dCurvilinear < scheme.Scheme 1 classdef Heat2dCurvilinear < scheme.Scheme
2 2
3 % Discretizes the Laplacian with variable coefficent, curvilinear, 3 % Discretizes the Laplacian with variable coefficent, curvilinear,
4 % in the Heat equation way (i.e., the discretization matrix is not necessarily 4 % in the Heat equation way (i.e., the discretization matrix is not necessarily
5 % symmetric) 5 % symmetric)
6 % u_t = div * (kappa * grad u ) 6 % u_t = div * (kappa * grad u )
7 % opSet should be cell array of opSets, one per dimension. This 7 % opSet should be cell array of opSets, one per dimension. This
8 % is useful if we have periodic BC in one direction. 8 % is useful if we have periodic BC in one direction.
9 9
10 properties 10 properties
11 m % Number of points in each direction, possibly a vector 11 m % Number of points in each direction, possibly a vector
27 27
28 H, Hi % Inner products 28 H, Hi % Inner products
29 e_l, e_r 29 e_l, e_r
30 d1_l, d1_r % Normal derivatives at the boundary 30 d1_l, d1_r % Normal derivatives at the boundary
31 alpha % Vector of borrowing constants 31 alpha % Vector of borrowing constants
32 32
33 % Boundary inner products 33 % Boundary inner products
34 H_boundary_l, H_boundary_r 34 H_boundary_l, H_boundary_r
35 35
36 % Metric coefficients 36 % Metric coefficients
37 b % Cell matrix of size dim x dim 37 b % Cell matrix of size dim x dim
38 J, Ji 38 J, Ji
39 beta % Cell array of scale factors 39 beta % Cell array of scale factors
107 xmax = max(ops{1}.x); 107 xmax = max(ops{1}.x);
108 ymax = max(ops{2}.x); 108 ymax = max(ops{2}.x);
109 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order); 109 opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order);
110 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order); 110 opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order);
111 D1Metric{1} = kron(opSetMetric{1}.D1, I{2}); 111 D1Metric{1} = kron(opSetMetric{1}.D1, I{2});
112 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); 112 D1Metric{2} = kron(I{1}, opSetMetric{2}.D1);
113 113
114 x_xi = D1Metric{1}*x; 114 x_xi = D1Metric{1}*x;
115 x_eta = D1Metric{2}*x; 115 x_eta = D1Metric{2}*x;
116 y_xi = D1Metric{1}*y; 116 y_xi = D1Metric{1}*y;
117 y_eta = D1Metric{2}*y; 117 y_eta = D1Metric{2}*y;
155 obj.d1_r{2} = kron(I{1},d1_r{2}); 155 obj.d1_r{2} = kron(I{1},d1_r{2});
156 156
157 % D2 coefficients 157 % D2 coefficients
158 kappa_coeff = cell(dim,dim); 158 kappa_coeff = cell(dim,dim);
159 for j = 1:dim 159 for j = 1:dim
160 obj.D2_kappa{j} = sparse(m_tot,m_tot); 160 obj.D2_kappa{j} = sparse(m_tot,m_tot);
161 kappa_coeff{j} = sparse(m_tot,1); 161 kappa_coeff{j} = sparse(m_tot,1);
162 for i = 1:dim 162 for i = 1:dim
163 kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa; 163 kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa;
164 end 164 end
165 end 165 end
268 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning) 268 function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning)
269 default_arg('type','Neumann'); 269 default_arg('type','Neumann');
270 default_arg('symmetric', false); 270 default_arg('symmetric', false);
271 default_arg('tuning',1.2); 271 default_arg('tuning',1.2);
272 272
273 % j is the coordinate direction of the boundary 273 % nj: outward unit normal component.
274 % nj: outward unit normal component.
275 % nj = -1 for west, south, bottom boundaries 274 % nj = -1 for west, south, bottom boundaries
276 % nj = 1 for east, north, top boundaries 275 % nj = 1 for east, north, top boundaries
277 [j, nj] = obj.get_boundary_number(boundary); 276 nj = obj.getBoundarySign(boundary);
278 switch nj 277
279 case 1 278 Hi = obj.Hi;
280 e = obj.e_r{j}; 279 [e, flux] = obj.getBoundaryOperator({'e', 'flux'}, boundary);
281 flux = obj.flux_r{j}; 280 H_gamma = obj.getBoundaryQuadrature(boundary);
282 H_gamma = obj.H_boundary_r{j}; 281 alpha = obj.getBoundaryBorrowing(boundary);
283 case -1
284 e = obj.e_l{j};
285 flux = obj.flux_l{j};
286 H_gamma = obj.H_boundary_l{j};
287 end
288 282
289 Hi = obj.Hi; 283 Hi = obj.Hi;
290 Ji = obj.Ji; 284 Ji = obj.Ji;
291 KAPPA = obj.KAPPA; 285 KAPPA = obj.KAPPA;
292 kappa_gamma = e'*KAPPA*e; 286 kappa_gamma = e'*KAPPA*e;
293 h = obj.h(j);
294 alpha = h*obj.alpha(j);
295 287
296 switch type 288 switch type
297 289
298 % Dirichlet boundary condition 290 % Dirichlet boundary condition
299 case {'D','d','dirichlet','Dirichlet'} 291 case {'D','d','dirichlet','Dirichlet'}
300 292
301 if ~symmetric 293 if ~symmetric
302 closure = -Ji*Hi*flux'*e*H_gamma*(e' ); 294 closure = -Ji*Hi*flux'*e*H_gamma*(e' );
303 penalty = Ji*Hi*flux'*e*H_gamma; 295 penalty = Ji*Hi*flux'*e*H_gamma;
304 else 296 else
305 closure = Ji*Hi*flux'*e*H_gamma*(e' )... 297 closure = Ji*Hi*flux'*e*H_gamma*(e' )...
306 -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ; 298 -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ;
307 penalty = -Ji*Hi*flux'*e*H_gamma ... 299 penalty = -Ji*Hi*flux'*e*H_gamma ...
308 +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma; 300 +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma;
309 end 301 end
310 302
311 % Normal flux boundary condition 303 % Normal flux boundary condition
312 case {'N','n','neumann','Neumann'} 304 case {'N','n','neumann','Neumann'}
313 closure = -Ji*Hi*e*H_gamma*(e'*flux ); 305 closure = -Ji*Hi*e*H_gamma*(e'*flux );
314 penalty = Ji*Hi*e*H_gamma; 306 penalty = Ji*Hi*e*H_gamma;
315 307
316 % Unknown boundary condition 308 % Unknown boundary condition
317 otherwise 309 otherwise
318 error('No such boundary condition: type = %s',type); 310 error('No such boundary condition: type = %s',type);
319 end 311 end
323 % u denotes the solution in the own domain 315 % u denotes the solution in the own domain
324 % v denotes the solution in the neighbour domain 316 % v denotes the solution in the neighbour domain
325 error('Interface not implemented'); 317 error('Interface not implemented');
326 end 318 end
327 319
328 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. 320 % Returns the boundary operator op for the boundary specified by the string boundary.
329 function [j, nj] = get_boundary_number(obj, boundary) 321 % op -- string or a cell array of strings
322 % boundary -- string
323 function varargout = getBoundaryOperator(obj, op, boundary)
324
325 if ~iscell(op)
326 op = {op};
327 end
328
329 for i = 1:numel(op)
330 switch op{i}
331 case 'e'
332 switch boundary
333 case 'w'
334 e = obj.e_l{1};
335 case 'e'
336 e = obj.e_r{1};
337 case 's'
338 e = obj.e_l{2};
339 case 'n'
340 e = obj.e_r{2};
341 otherwise
342 error('No such boundary: boundary = %s',boundary);
343 end
344 varargout{i} = e;
345
346 case 'd'
347 switch boundary
348 case 'w'
349 d = obj.d1_l{1};
350 case 'e'
351 d = obj.d1_r{1};
352 case 's'
353 d = obj.d1_l{2};
354 case 'n'
355 d = obj.d1_r{2};
356 otherwise
357 error('No such boundary: boundary = %s',boundary);
358 end
359 varargout{i} = d;
360
361 case 'flux'
362 switch boundary
363 case 'w'
364 flux = obj.flux_l{1};
365 case 'e'
366 flux = obj.flux_r{1};
367 case 's'
368 flux = obj.flux_l{2};
369 case 'n'
370 flux = obj.flux_r{2};
371 otherwise
372 error('No such boundary: boundary = %s',boundary);
373 end
374 varargout{i} = flux;
375 end
376 end
377 end
378
379 % Returns square boundary quadrature matrix, of dimension
380 % corresponding to the number of boundary points
381 %
382 % boundary -- string
383 function H_b = getBoundaryQuadrature(obj, boundary)
330 384
331 switch boundary 385 switch boundary
332 case {'w','W','west','West', 'e', 'E', 'east', 'East'} 386 case 'w'
333 j = 1; 387 H_b = obj.H_boundary_l{1};
334 case {'s','S','south','South', 'n', 'N', 'north', 'North'} 388 case 'e'
335 j = 2; 389 H_b = obj.H_boundary_r{1};
390 case 's'
391 H_b = obj.H_boundary_l{2};
392 case 'n'
393 H_b = obj.H_boundary_r{2};
336 otherwise 394 otherwise
337 error('No such boundary: boundary = %s',boundary); 395 error('No such boundary: boundary = %s',boundary);
338 end 396 end
339 397 end
398
399 % Returns the boundary sign. The right boundary is considered the positive boundary
400 % boundary -- string
401 function s = getBoundarySign(obj, boundary)
340 switch boundary 402 switch boundary
341 case {'w','W','west','West','s','S','south','South'} 403 case {'e','n'}
342 nj = -1; 404 s = 1;
343 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 405 case {'w','s'}
344 nj = 1; 406 s = -1;
345 end
346 end
347
348 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
349 function [return_op] = get_boundary_operator(obj, op, boundary)
350
351 switch boundary
352 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
353 j = 1;
354 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
355 j = 2;
356 otherwise 407 otherwise
357 error('No such boundary: boundary = %s',boundary); 408 error('No such boundary: boundary = %s',boundary);
358 end 409 end
359 410 end
360 switch op 411
361 case 'e' 412 % Returns borrowing constant gamma*h
362 switch boundary 413 % boundary -- string
363 case {'w','W','west','West','s','S','south','South'} 414 function gamm = getBoundaryBorrowing(obj, boundary)
364 return_op = obj.e_l{j}; 415 switch boundary
365 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} 416 case {'w','e'}
366 return_op = obj.e_r{j}; 417 gamm = obj.h(1)*obj.alpha(1);
367 end 418 case {'s','n'}
368 case 'd' 419 gamm = obj.h(2)*obj.alpha(2);
369 switch boundary
370 case {'w','W','west','West','s','S','south','South'}
371 return_op = obj.d1_l{j};
372 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
373 return_op = obj.d1_r{j};
374 end
375 otherwise 420 otherwise
376 error(['No such operator: operatr = ' op]); 421 error('No such boundary: boundary = %s',boundary);
377 end 422 end
378
379 end 423 end
380 424
381 function N = size(obj) 425 function N = size(obj)
382 N = prod(obj.m); 426 N = prod(obj.m);
383 end 427 end