Mercurial > repos > public > sbplib

comparison +scheme/LaplaceCurvilinear.m @ 938:97291e1bd57c feature/utux2D

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Make LaplCurv send only number of grid points to InterpOpsXX. Remove coordinates from get_boundary_ops method.
author Martin Almquist <malmquist@stanford.edu>
date Tue, 04 Dec 2018 14:07:38 -0800
parents a205264664cd
children 31186559236d
comparison
equal deleted inserted replaced
937:ed8c98c4d479 938:97291e1bd57c
334 tuning = type.tuning; 334 tuning = type.tuning;
335 335
336 336
337 % u denotes the solution in the own domain 337 % u denotes the solution in the own domain
338 % v denotes the solution in the neighbour domain 338 % v denotes the solution in the neighbour domain
339 [e_u, d_u, gamm_u, H_b_u, I_u, ~, X_u] = obj.get_boundary_ops(boundary); 339 [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary);
340 [e_v, d_v, gamm_v, H_b_v, I_v, ~, X_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 340 [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
341 341
342 % Extract the coordinate that varies along the interface 342 % Find the number of grid points along the interface
343 switch boundary 343 m_u = size(e_u, 2);
344 case {'e','w'} 344 m_v = size(e_v, 2);
345 x_u = X_u(:,2);
346 case {'s', 'n'}
347 x_u = X_u(:,1);
348 end
349
350 switch neighbour_boundary
351 case {'e','w'}
352 x_v = X_v(:,2);
353 case {'s', 'n'}
354 x_v = X_v(:,1);
355 end
356 345
357 Hi = obj.Hi; 346 Hi = obj.Hi;
358 a = obj.a; 347 a = obj.a;
359 348
360 u = obj; 349 u = obj;
371 tau_u = tuning * spdiag(tau_u); 360 tau_u = tuning * spdiag(tau_u);
372 tau_v = tuning * spdiag(tau_v); 361 tau_v = tuning * spdiag(tau_v);
373 beta_u = tau_v; 362 beta_u = tau_v;
374 363
375 % Build interpolation operators 364 % Build interpolation operators
376 intOps = interpOpSet(x_u, x_v, obj.order, neighbour_scheme.order); 365 intOps = interpOpSet(m_u, m_v, obj.order, neighbour_scheme.order);
377 Iu2v = intOps.Iu2v; 366 Iu2v = intOps.Iu2v;
378 Iv2u = intOps.Iv2u; 367 Iv2u = intOps.Iv2u;
379 368
380 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ... 369 closure = a*Hi*e_u*tau_u*H_b_u*e_u' + ...
381 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ... 370 a*Hi*e_u*H_b_u*Iv2u.bad*beta_u*Iu2v.good*e_u' + ...
391 380
392 % Returns the boundary ops and sign for the boundary specified by the string boundary. 381 % Returns the boundary ops and sign for the boundary specified by the string boundary.
393 % The right boundary is considered the positive boundary 382 % The right boundary is considered the positive boundary
394 % 383 %
395 % I -- the indices of the boundary points in the grid matrix 384 % I -- the indices of the boundary points in the grid matrix
396 % X -- coordinates along the boundary; 385 function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary)
397 function [e, d, gamm, H_b, I, X, X_logic] = get_boundary_ops(obj, boundary)
398 386
399 % gridMatrix = zeros(obj.m(2),obj.m(1)); 387 % gridMatrix = zeros(obj.m(2),obj.m(1));
400 % gridMatrix(:) = 1:numel(gridMatrix); 388 % gridMatrix(:) = 1:numel(gridMatrix);
401
402 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); 389 ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
403 390
404 switch boundary 391 switch boundary
405 case 'w' 392 case 'w'
406 e = obj.e_w; 393 e = obj.e_w;
424 I = ind(:,end)'; 411 I = ind(:,end)';
425 otherwise 412 otherwise
426 error('No such boundary: boundary = %s',boundary); 413 error('No such boundary: boundary = %s',boundary);
427 end 414 end
428 415
429 X = obj.grid.getBoundary(boundary);
430 if isa(obj.grid, 'grid.Curvilinear')
431 X_logic = obj.grid.logic.getBoundary(boundary);
432 else
433 % Cartesian physical coordinates are also logical coordinates
434 X_logic = X;
435 end
436
437 switch boundary 416 switch boundary
438 case {'w','e'} 417 case {'w','e'}
439 gamm = obj.gamm_u; 418 gamm = obj.gamm_u;
440 case {'s','n'} 419 case {'s','n'}
441 gamm = obj.gamm_v; 420 gamm = obj.gamm_v;