Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 1197:433c89bf19e0 feature/rv

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Wed, 07 Aug 2019 15:23:42 +0200
parents 84200bbae101
children 5ec23b9bf360 d1dad4fbfe22
comparison
equal deleted inserted replaced
1196:f6c571d8f22f 1197:433c89bf19e0
236 % neighbour_boundary is a string specifying which boundary to interface to. 236 % neighbour_boundary is a string specifying which boundary to interface to.
237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) 237 function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
238 default_arg('type','neumann'); 238 default_arg('type','neumann');
239 default_arg('parameter', []); 239 default_arg('parameter', []);
240 240
241 [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); 241 e = obj.getBoundaryOperator('e', boundary);
242 d = obj.getBoundaryOperator('d', boundary);
243 H_b = obj.getBoundaryQuadrature(boundary);
244 gamm = obj.getBoundaryBorrowing(boundary);
242 switch type 245 switch type
243 % Dirichlet boundary condition 246 % Dirichlet boundary condition
244 case {'D','d','dirichlet'} 247 case {'D','d','dirichlet'}
245 tuning = 1.2; 248 tuning = 1.2;
246 % tuning = 20.2; 249 % tuning = 20.2;
296 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) 299 function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type)
297 tuning = type.tuning; 300 tuning = type.tuning;
298 301
299 % u denotes the solution in the own domain 302 % u denotes the solution in the own domain
300 % v denotes the solution in the neighbour domain 303 % v denotes the solution in the neighbour domain
301 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 304 e_u = obj.getBoundaryOperator('e', boundary);
302 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 305 d_u = obj.getBoundaryOperator('d', boundary);
306 H_b_u = obj.getBoundaryQuadrature(boundary);
307 I_u = obj.getBoundaryIndices(boundary);
308 gamm_u = obj.getBoundaryBorrowing(boundary);
309
310 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
311 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
312 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
313 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
314 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
303 315
304 u = obj; 316 u = obj;
305 v = neighbour_scheme; 317 v = neighbour_scheme;
306 318
307 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; 319 b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2;
334 tuning = type.tuning; 346 tuning = type.tuning;
335 347
336 348
337 % u denotes the solution in the own domain 349 % u denotes the solution in the own domain
338 % v denotes the solution in the neighbour domain 350 % v denotes the solution in the neighbour domain
339 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); 351 e_u = obj.getBoundaryOperator('e', boundary);
340 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 352 d_u = obj.getBoundaryOperator('d', boundary);
353 H_b_u = obj.getBoundaryQuadrature(boundary);
354 I_u = obj.getBoundaryIndices(boundary);
355 gamm_u = obj.getBoundaryBorrowing(boundary);
356
357 e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
358 d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
359 H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
360 I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
361 gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
362
341 363
342 % Find the number of grid points along the interface 364 % Find the number of grid points along the interface
343 m_u = size(e_u, 2); 365 m_u = size(e_u, 2);
344 m_v = size(e_v, 2); 366 m_v = size(e_v, 2);
345 367
376 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ... 398 -a*1/2*Hi*d_u*H_b_u*Iv2u.good*e_v' + ...
377 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v'; 399 -a*1/2*Hi*e_u*H_b_u*Iv2u.bad*d_v';
378 400
379 end 401 end
380 402
381 % Returns the boundary ops and sign for the boundary specified by the string boundary. 403 % Returns the boundary operator op for the boundary specified by the string boundary.
382 % The right boundary is considered the positive boundary 404 % op -- string
405 % boundary -- string
406 function o = getBoundaryOperator(obj, op, boundary)
407 assertIsMember(op, {'e', 'd'})
408 assertIsMember(boundary, {'w', 'e', 's', 'n'})
409
410 o = obj.([op, '_', boundary]);
411 end
412
413 % Returns square boundary quadrature matrix, of dimension
414 % corresponding to the number of boundary points
383 % 415 %
384 % I -- the indices of the boundary points in the grid matrix 416 % boundary -- string
385 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) 417 function H_b = getBoundaryQuadrature(obj, boundary)
418 assertIsMember(boundary, {'w', 'e', 's', 'n'})
419
420 H_b = obj.(['H_', boundary]);
421 end
422
423 % Returns the indices of the boundary points in the grid matrix
424 % boundary -- string
425 function I = getBoundaryIndices(obj, boundary)
426 assertIsMember(boundary, {'w', 'e', 's', 'n'})
427
386 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 428 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
387
388 switch boundary 429 switch boundary
389 case 'w' 430 case 'w'
390 e = obj.e_w;
391 d = obj.d_w;
392 H_b = obj.H_w;
393 I = ind(1,:); 431 I = ind(1,:);
394 case 'e' 432 case 'e'
395 e = obj.e_e;
396 d = obj.d_e;
397 H_b = obj.H_e;
398 I = ind(end,:); 433 I = ind(end,:);
399 case 's' 434 case 's'
400 e = obj.e_s;
401 d = obj.d_s;
402 H_b = obj.H_s;
403 I = ind(:,1)'; 435 I = ind(:,1)';
404 case 'n' 436 case 'n'
405 e = obj.e_n;
406 d = obj.d_n;
407 H_b = obj.H_n;
408 I = ind(:,end)'; 437 I = ind(:,end)';
409 otherwise 438 end
410 error('No such boundary: boundary = %s',boundary); 439 end
411 end 440
441 % Returns borrowing constant gamma
442 % boundary -- string
443 function gamm = getBoundaryBorrowing(obj, boundary)
444 assertIsMember(boundary, {'w', 'e', 's', 'n'})
412 445
413 switch boundary 446 switch boundary
414 case {'w','e'} 447 case {'w','e'}
415 gamm = obj.gamm_u; 448 gamm = obj.gamm_u;
416 case {'s','n'} 449 case {'s','n'}