Mercurial > repos > public > sbplib
changeset 997:78db023a7fe3 feature/getBoundaryOp
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 | a0b3161e44f3 |
| children | 2b1b944deae1 |
| files | +scheme/Beam2d.m +scheme/Heat2dCurvilinear.m +scheme/Heat2dVariable.m +scheme/Hypsyst2d.m +scheme/Hypsyst2dCurve.m +scheme/LaplaceCurvilinear.m +scheme/Schrodinger2d.m +scheme/Utux2d.m +scheme/Wave2d.m |
| diffstat | 9 files changed, 683 insertions(+), 246 deletions(-) [+] |
line wrap: on
line diff
--- a/+scheme/Beam2d.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Beam2d.m Sat Jan 12 11:57:50 2019 -0800 @@ -121,7 +121,10 @@ default_arg('type','dn'); default_arg('data',0); - [e,d1,d2,d3,s,gamm,delt,halfnorm_inv] = obj.get_boundary_ops(boundary); + [e, d1, d2, d3] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary); + s = obj.getBoundarySign(boundary); + [gamm, delt] = obj.getBoundaryBorrowing(boundary); + halfnorm_inv = obj.getHalfnormInv(boundary); switch type % Dirichlet-neumann boundary condition @@ -164,8 +167,14 @@ function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary, type) % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - [e_u,d1_u,d2_u,d3_u,s_u,gamm_u,delt_u, halfnorm_inv] = obj.get_boundary_ops(boundary); - [e_v,d1_v,d2_v,d3_v,s_v,gamm_v,delt_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u, d1_u, d2_u, d3_u] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary); + s_u = obj.getBoundarySign(boundary); + [gamm_u, delt_u] = obj.getBoundaryBorrowing(boundary); + halfnorm_inv = obj.getHalfnormInv(boundary); + + [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary); + s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); + [gamm_v, delt_v] = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); tuning = 2; @@ -192,46 +201,141 @@ penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); end - % Ruturns the boundary ops and sign for the boundary specified by the string boundary. - % The right boundary is considered the positive boundary - function [e,d1,d2,d3,s,gamm, delt, halfnorm_inv] = get_boundary_ops(obj,boundary) + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_w; + case 'e' + e = obj.e_e; + case 's' + e = obj.e_s; + case 'n' + e = obj.e_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd1' + switch boundary + case 'w' + d1 = obj.d1_w; + case 'e' + d1 = obj.d1_e; + case 's' + d1 = obj.d1_s; + case 'n' + d1 = obj.d1_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d1; + end + + case 'd2' + switch boundary + case 'w' + d2 = obj.d2_w; + case 'e' + d2 = obj.d2_e; + case 's' + d2 = obj.d2_s; + case 'n' + d2 = obj.d2_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d2; + end + + case 'd3' + switch boundary + case 'w' + d3 = obj.d3_w; + case 'e' + d3 = obj.d3_e; + case 's' + d3 = obj.d3_s; + case 'n' + d3 = obj.d3_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d3; + end + end + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) + switch boundary case 'w' - e = obj.e_w; - d1 = obj.d1_w; - d2 = obj.d2_w; - d3 = obj.d3_w; + H_b = obj.H_y; + case 'e' + H_b = obj.H_y; + case 's' + H_b = obj.H_x; + case 'n' + H_b = obj.H_x; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + switch boundary + case {'e','n'} + s = 1; + case {'w','s'} s = -1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns the halfnorm_inv used in SATs. TODO: better notation + function Hinv = getHalfnormInv(obj, boundary) + switch boundary + case 'w' + Hinv = obj.Hix; + case 'e' + Hinv = obj.Hix; + case 's' + Hinv = obj.Hiy; + case 'n' + Hinv = obj.Hiy; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns borrowing constant gamma + % boundary -- string + function [gamm, delt] = getBoundaryBorrowing(obj, boundary) + switch boundary + case {'w','e'} gamm = obj.gamm_x; delt = obj.delt_x; - halfnorm_inv = obj.Hix; - case 'e' - e = obj.e_e; - d1 = obj.d1_e; - d2 = obj.d2_e; - d3 = obj.d3_e; - s = 1; - gamm = obj.gamm_x; - delt = obj.delt_x; - halfnorm_inv = obj.Hix; - case 's' - e = obj.e_s; - d1 = obj.d1_s; - d2 = obj.d2_s; - d3 = obj.d3_s; - s = -1; + case {'s','n'} gamm = obj.gamm_y; delt = obj.delt_y; - halfnorm_inv = obj.Hiy; - case 'n' - e = obj.e_n; - d1 = obj.d1_n; - d2 = obj.d2_n; - d3 = obj.d3_n; - s = 1; - gamm = obj.gamm_y; - delt = obj.delt_y; - halfnorm_inv = obj.Hiy; otherwise error('No such boundary: boundary = %s',boundary); end
--- a/+scheme/Heat2dCurvilinear.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Heat2dCurvilinear.m Sat Jan 12 11:57:50 2019 -0800 @@ -1,9 +1,9 @@ classdef Heat2dCurvilinear < scheme.Scheme % Discretizes the Laplacian with variable coefficent, curvilinear, -% in the Heat equation way (i.e., the discretization matrix is not necessarily +% in the Heat equation way (i.e., the discretization matrix is not necessarily % symmetric) -% u_t = div * (kappa * grad u ) +% u_t = div * (kappa * grad u ) % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -29,9 +29,9 @@ e_l, e_r d1_l, d1_r % Normal derivatives at the boundary alpha % Vector of borrowing constants - + % Boundary inner products - H_boundary_l, H_boundary_r + H_boundary_l, H_boundary_r % Metric coefficients b % Cell matrix of size dim x dim @@ -109,7 +109,7 @@ opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order); opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order); D1Metric{1} = kron(opSetMetric{1}.D1, I{2}); - D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); + D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); x_xi = D1Metric{1}*x; x_eta = D1Metric{2}*x; @@ -157,7 +157,7 @@ % D2 coefficients kappa_coeff = cell(dim,dim); for j = 1:dim - obj.D2_kappa{j} = sparse(m_tot,m_tot); + obj.D2_kappa{j} = sparse(m_tot,m_tot); kappa_coeff{j} = sparse(m_tot,1); for i = 1:dim kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa; @@ -270,28 +270,20 @@ default_arg('symmetric', false); default_arg('tuning',1.2); - % j is the coordinate direction of the boundary - % nj: outward unit normal component. + % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries % nj = 1 for east, north, top boundaries - [j, nj] = obj.get_boundary_number(boundary); - switch nj - case 1 - e = obj.e_r{j}; - flux = obj.flux_r{j}; - H_gamma = obj.H_boundary_r{j}; - case -1 - e = obj.e_l{j}; - flux = obj.flux_l{j}; - H_gamma = obj.H_boundary_l{j}; - end + nj = obj.getBoundarySign(boundary); + + Hi = obj.Hi; + [e, flux] = obj.getBoundaryOperator({'e', 'flux'}, boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); + alpha = obj.getBoundaryBorrowing(boundary); Hi = obj.Hi; Ji = obj.Ji; KAPPA = obj.KAPPA; - kappa_gamma = e'*KAPPA*e; - h = obj.h(j); - alpha = h*obj.alpha(j); + kappa_gamma = e'*KAPPA*e; switch type @@ -299,19 +291,19 @@ case {'D','d','dirichlet','Dirichlet'} if ~symmetric - closure = -Ji*Hi*flux'*e*H_gamma*(e' ); + closure = -Ji*Hi*flux'*e*H_gamma*(e' ); penalty = Ji*Hi*flux'*e*H_gamma; else closure = Ji*Hi*flux'*e*H_gamma*(e' )... - -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ; + -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ; penalty = -Ji*Hi*flux'*e*H_gamma ... +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma; end % Normal flux boundary condition case {'N','n','neumann','Neumann'} - closure = -Ji*Hi*e*H_gamma*(e'*flux ); - penalty = Ji*Hi*e*H_gamma; + closure = -Ji*Hi*e*H_gamma*(e'*flux ); + penalty = Ji*Hi*e*H_gamma; % Unknown boundary condition otherwise @@ -325,57 +317,109 @@ error('Interface not implemented'); end - % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. - function [j, nj] = get_boundary_number(obj, boundary) + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) - switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} - j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} - j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); + if ~iscell(op) + op = {op}; end - switch boundary - case {'w','W','west','West','s','S','south','South'} - nj = -1; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - nj = 1; + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_l{1}; + case 'e' + e = obj.e_r{1}; + case 's' + e = obj.e_l{2}; + case 'n' + e = obj.e_r{2}; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd' + switch boundary + case 'w' + d = obj.d1_l{1}; + case 'e' + d = obj.d1_r{1}; + case 's' + d = obj.d1_l{2}; + case 'n' + d = obj.d1_r{2}; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d; + + case 'flux' + switch boundary + case 'w' + flux = obj.flux_l{1}; + case 'e' + flux = obj.flux_r{1}; + case 's' + flux = obj.flux_l{2}; + case 'n' + flux = obj.flux_r{2}; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = flux; + end end end - % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. - function [return_op] = get_boundary_operator(obj, op, boundary) + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} - j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} - j = 2; + case 'w' + H_b = obj.H_boundary_l{1}; + case 'e' + H_b = obj.H_boundary_r{1}; + case 's' + H_b = obj.H_boundary_l{2}; + case 'n' + H_b = obj.H_boundary_r{2}; otherwise error('No such boundary: boundary = %s',boundary); end + end - switch op - case 'e' - switch boundary - case {'w','W','west','West','s','S','south','South'} - return_op = obj.e_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - return_op = obj.e_r{j}; - end - case 'd' - switch boundary - case {'w','W','west','West','s','S','south','South'} - return_op = obj.d1_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - return_op = obj.d1_r{j}; - end + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + switch boundary + case {'e','n'} + s = 1; + case {'w','s'} + s = -1; otherwise - error(['No such operator: operatr = ' op]); + error('No such boundary: boundary = %s',boundary); end + end + % Returns borrowing constant gamma*h + % boundary -- string + function gamm = getBoundaryBorrowing(obj, boundary) + switch boundary + case {'w','e'} + gamm = obj.h(1)*obj.alpha(1); + case {'s','n'} + gamm = obj.h(2)*obj.alpha(2); + otherwise + error('No such boundary: boundary = %s',boundary); + end end function N = size(obj)
--- a/+scheme/Heat2dVariable.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Heat2dVariable.m Sat Jan 12 11:57:50 2019 -0800 @@ -1,9 +1,9 @@ classdef Heat2dVariable < scheme.Scheme % Discretizes the Laplacian with variable coefficent, -% In the Heat equation way (i.e., the discretization matrix is not necessarily +% In the Heat equation way (i.e., the discretization matrix is not necessarily % symmetric) -% u_t = div * (kappa * grad u ) +% u_t = div * (kappa * grad u ) % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -29,7 +29,7 @@ e_l, e_r d1_l, d1_r % Normal derivatives at the boundary alpha % Vector of borrowing constants - + H_boundary % Boundary inner products end @@ -162,26 +162,18 @@ default_arg('symmetric', false); default_arg('tuning',1.2); - % j is the coordinate direction of the boundary - % nj: outward unit normal component. + % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries % nj = 1 for east, north, top boundaries - [j, nj] = obj.get_boundary_number(boundary); - switch nj - case 1 - e = obj.e_r; - d = obj.d1_r; - case -1 - e = obj.e_l; - d = obj.d1_l; - end + nj = obj.getBoundarySign(boundary); Hi = obj.Hi; - H_gamma = obj.H_boundary{j}; + [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); + alpha = obj.getBoundaryBorrowing(boundary); + KAPPA = obj.KAPPA; - kappa_gamma = e{j}'*KAPPA*e{j}; - h = obj.h(j); - alpha = h*obj.alpha(j); + kappa_gamma = e'*KAPPA*e; switch type @@ -189,19 +181,19 @@ case {'D','d','dirichlet','Dirichlet'} if ~symmetric - closure = -nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' ); - penalty = nj*Hi*d{j}*kappa_gamma*H_gamma; + closure = -nj*Hi*d*kappa_gamma*H_gamma*(e' ); + penalty = nj*Hi*d*kappa_gamma*H_gamma; else - closure = nj*Hi*d{j}*kappa_gamma*H_gamma*(e{j}' )... - -tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma*(e{j}' ) ; - penalty = -nj*Hi*d{j}*kappa_gamma*H_gamma ... - +tuning*2/alpha*Hi*e{j}*kappa_gamma*H_gamma; + closure = nj*Hi*d*kappa_gamma*H_gamma*(e' )... + -tuning*2/alpha*Hi*e*kappa_gamma*H_gamma*(e' ) ; + penalty = -nj*Hi*d*kappa_gamma*H_gamma ... + +tuning*2/alpha*Hi*e*kappa_gamma*H_gamma; end % Free boundary condition case {'N','n','neumann','Neumann'} - closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); - penalty = Hi*e{j}*kappa_gamma*H_gamma; + closure = -nj*Hi*e*kappa_gamma*H_gamma*(d' ); + penalty = Hi*e*kappa_gamma*H_gamma; % penalty is for normal derivative and not for derivative, hence the sign. % Unknown boundary condition @@ -216,57 +208,94 @@ error('Interface not implemented'); end - % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. - function [j, nj] = get_boundary_number(obj, boundary) + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) - switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} - j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} - j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); + if ~iscell(op) + op = {op}; end - switch boundary - case {'w','W','west','West','s','S','south','South'} - nj = -1; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - nj = 1; + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_l{1}; + case 'e' + e = obj.e_r{1}; + case 's' + e = obj.e_l{2}; + case 'n' + e = obj.e_r{2}; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd' + switch boundary + case 'w' + d = obj.d1_l{1}; + case 'e' + d = obj.d1_r{1}; + case 's' + d = obj.d1_l{2}; + case 'n' + d = obj.d1_r{2}; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d; + end end end - % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. - function [return_op] = get_boundary_operator(obj, op, boundary) + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} - j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} - j = 2; + case 'w' + H_b = obj.H_boundary{1}; + case 'e' + H_b = obj.H_boundary{1}; + case 's' + H_b = obj.H_boundary{2}; + case 'n' + H_b = obj.H_boundary{2}; otherwise error('No such boundary: boundary = %s',boundary); end + end - switch op - case 'e' - switch boundary - case {'w','W','west','West','s','S','south','South'} - return_op = obj.e_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - return_op = obj.e_r{j}; - end - case 'd' - switch boundary - case {'w','W','west','West','s','S','south','South'} - return_op = obj.d1_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - return_op = obj.d1_r{j}; - end + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + switch boundary + case {'e','n'} + s = 1; + case {'w','s'} + s = -1; otherwise - error(['No such operator: operatr = ' op]); + error('No such boundary: boundary = %s',boundary); end + end + % Returns borrowing constant gamma*h + % boundary -- string + function gamm = getBoundaryBorrowing(obj, boundary) + switch boundary + case {'w','e'} + gamm = obj.h(1)*obj.alpha(1); + case {'s','n'} + gamm = obj.h(2)*obj.alpha(2); + otherwise + error('No such boundary: boundary = %s',boundary); + end end function N = size(obj)
--- a/+scheme/Hypsyst2d.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Hypsyst2d.m Sat Jan 12 11:57:50 2019 -0800 @@ -186,10 +186,10 @@ params = obj.params; x = obj.x; y = obj.y; + e_ = obj.getBoundaryOperator('e', boundary); switch boundary case {'w','W','west'} - e_ = obj.e_w; mat = obj.A; boundPos = 'l'; Hi = obj.Hxi; @@ -197,7 +197,6 @@ L = obj.evaluateCoefficientMatrix(L,x(1),y); side = max(length(y)); case {'e','E','east'} - e_ = obj.e_e; mat = obj.A; boundPos = 'r'; Hi = obj.Hxi; @@ -205,7 +204,6 @@ L = obj.evaluateCoefficientMatrix(L,x(end),y); side = max(length(y)); case {'s','S','south'} - e_ = obj.e_s; mat = obj.B; boundPos = 'l'; Hi = obj.Hyi; @@ -213,7 +211,6 @@ L = obj.evaluateCoefficientMatrix(L,x,y(1)); side = max(length(x)); case {'n','N','north'} - e_ = obj.e_n; mat = obj.B; boundPos = 'r'; Hi = obj.Hyi; @@ -297,5 +294,56 @@ signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; end + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_w; + case 'e' + e = obj.e_e; + case 's' + e = obj.e_s; + case 'n' + e = obj.e_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + end + end + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) + + e = obj.getBoundaryOperator('e', boundary); + + switch boundary + case 'w' + H_b = inv(e'*obj.Hyi*e); + case 'e' + H_b = inv(e'*obj.Hyi*e); + case 's' + H_b = inv(e'*obj.Hxi*e); + case 'n' + H_b = inv(e'*obj.Hxi*e); + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + end end \ No newline at end of file
--- a/+scheme/Hypsyst2dCurve.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Hypsyst2dCurve.m Sat Jan 12 11:57:50 2019 -0800 @@ -169,31 +169,28 @@ Y = obj.Y; xi = obj.xi; eta = obj.eta; + e_ = obj.getBoundaryOperator('e', boundary); switch boundary case {'w','W','west'} - e_ = obj.e_w; mat = obj.Ahat; boundPos = 'l'; Hi = obj.Hxii; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w)); side = max(length(eta)); case {'e','E','east'} - e_ = obj.e_e; mat = obj.Ahat; boundPos = 'r'; Hi = obj.Hxii; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e)); side = max(length(eta)); case {'s','S','south'} - e_ = obj.e_s; mat = obj.Bhat; boundPos = 'l'; Hi = obj.Hetai; [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s)); side = max(length(xi)); case {'n','N','north'} - e_ = obj.e_n; mat = obj.Bhat; boundPos = 'r'; Hi = obj.Hetai; @@ -374,5 +371,58 @@ Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)]; signVec = [sum(poseig),sum(zeroeig),sum(negeig)]; end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_w; + case 'e' + e = obj.e_e; + case 's' + e = obj.e_s; + case 'n' + e = obj.e_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + end + end + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) + + e = obj.getBoundaryOperator('e', boundary); + + switch boundary + case 'w' + H_b = inv(e'*obj.Hetai*e); + case 'e' + H_b = inv(e'*obj.Hetai*e); + case 's' + H_b = inv(e'*obj.Hxii*e); + case 'n' + H_b = inv(e'*obj.Hxii*e); + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + end end \ No newline at end of file
--- a/+scheme/LaplaceCurvilinear.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/LaplaceCurvilinear.m Sat Jan 12 11:57:50 2019 -0800 @@ -437,7 +437,6 @@ varargout{i} = d; end end - end % Returns square boundary quadrature matrix, of dimension
--- a/+scheme/Schrodinger2d.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Schrodinger2d.m Sat Jan 12 11:57:50 2019 -0800 @@ -162,35 +162,26 @@ default_arg('type','Neumann'); default_arg('parameter', []); - % j is the coordinate direction of the boundary % nj: outward unit normal component. % nj = -1 for west, south, bottom boundaries % nj = 1 for east, north, top boundaries - [j, nj] = obj.get_boundary_number(boundary); - switch nj - case 1 - e = obj.e_r; - d = obj.d1_r; - case -1 - e = obj.e_l; - d = obj.d1_l; - end - + nj = obj.getBoundarySign(boundary); + [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); Hi = obj.Hi; - H_gamma = obj.H_boundary{j}; - a = e{j}'*obj.a*e{j}; + a = e'*obj.a*e; switch type % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - closure = nj*Hi*d{j}*a*1i*H_gamma*(e{j}' ); - penalty = -nj*Hi*d{j}*a*1i*H_gamma; + closure = nj*Hi*d*a*1i*H_gamma*(e' ); + penalty = -nj*Hi*d*a*1i*H_gamma; % Free boundary condition case {'N','n','neumann','Neumann'} - closure = -nj*Hi*e{j}*a*1i*H_gamma*(d{j}' ); - penalty = nj*Hi*e{j}*a*1i*H_gamma; + closure = -nj*Hi*e*a*1i*H_gamma*(d' ); + penalty = nj*Hi*e*a*1i*H_gamma; % Unknown boundary condition otherwise @@ -221,13 +212,14 @@ % v denotes the solution in the neighbour domain % Get boundary operators - [e_neighbour, d_neighbour] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - [e, d, H_gamma] = obj.get_boundary_ops(boundary); + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); Hi = obj.Hi; a = obj.a; % Get outward unit normal component - [~, n] = obj.get_boundary_number(boundary); + n = obj.getBoundarySign(boundary); Hi = obj.Hi; sigma = -n*1i*a/2; @@ -247,13 +239,14 @@ % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - [e_v, d_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - [e_u, d_u, H_gamma] = obj.get_boundary_ops(boundary); + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); Hi = obj.Hi; a = obj.a; % Get outward unit normal component - [~, n] = obj.get_boundary_number(boundary); + n = obj.getBoundarySign(boundary); % Find the number of grid points along the interface m_u = size(e_u, 2); @@ -293,32 +286,83 @@ end end - % Returns the boundary ops and sign for the boundary specified by the string boundary. - % The right boundary is considered the positive boundary - function [e, d, H_b] = get_boundary_ops(obj, boundary) + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_w; + case 'e' + e = obj.e_e; + case 's' + e = obj.e_s; + case 'n' + e = obj.e_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd' + switch boundary + case 'w' + d = obj.d_w; + case 'e' + d = obj.d_e; + case 's' + d = obj.d_s; + case 'n' + d = obj.d_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d; + end + end + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) switch boundary case 'w' - e = obj.e_w; - d = obj.d_w; H_b = obj.H_boundary{1}; case 'e' - e = obj.e_e; - d = obj.d_e; H_b = obj.H_boundary{1}; case 's' - e = obj.e_s; - d = obj.d_s; H_b = obj.H_boundary{2}; case 'n' - e = obj.e_n; - d = obj.d_n; H_b = obj.H_boundary{2}; otherwise error('No such boundary: boundary = %s',boundary); end end + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + switch boundary + case {'e','n'} + s = 1; + case {'w','s'} + s = -1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + function N = size(obj) N = prod(obj.m); end
--- a/+scheme/Utux2d.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Utux2d.m Sat Jan 12 11:57:50 2019 -0800 @@ -139,16 +139,7 @@ couplingType = type.couplingType; % Get neighbour boundary operator - switch neighbour_boundary - case {'e','E','east','East'} - e_neighbour = neighbour_scheme.e_e; - case {'w','W','west','West'} - e_neighbour = neighbour_scheme.e_w; - case {'n','N','north','North'} - e_neighbour = neighbour_scheme.e_n; - case {'s','S','south','South'} - e_neighbour = neighbour_scheme.e_s; - end + e_neighbour = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); switch couplingType @@ -197,16 +188,7 @@ interpolationDamping = type.interpolationDamping; % Get neighbour boundary operator - switch neighbour_boundary - case {'e','E','east','East'} - e_neighbour = neighbour_scheme.e_e; - case {'w','W','west','West'} - e_neighbour = neighbour_scheme.e_w; - case {'n','N','north','North'} - e_neighbour = neighbour_scheme.e_n; - case {'s','S','south','South'} - e_neighbour = neighbour_scheme.e_s; - end + e_neighbour = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); switch couplingType @@ -290,6 +272,56 @@ end + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_w; + case 'e' + e = obj.e_e; + case 's' + e = obj.e_s; + case 'n' + e = obj.e_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + end + end + + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) + + switch boundary + case 'w' + H_b = obj.H_y; + case 'e' + H_b = obj.H_y; + case 's' + H_b = obj.H_x; + case 'n' + H_b = obj.H_x; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + function N = size(obj) N = obj.m; end
--- a/+scheme/Wave2d.m Tue Jan 08 11:51:24 2019 +0100 +++ b/+scheme/Wave2d.m Sat Jan 12 11:57:50 2019 -0800 @@ -106,7 +106,10 @@ default_arg('type','neumann'); default_arg('data',0); - [e,d,s,gamm,halfnorm_inv] = obj.get_boundary_ops(boundary); + [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + gamm = obj.getBoundaryBorrowing(boundary); + s = obj.getBoundarySign(boundary); + halfnorm_inv = obj.getHalfnormInv(boundary); switch type % Dirichlet boundary condition @@ -164,6 +167,15 @@ [e_u,d_u,s_u,gamm_u, halfnorm_inv] = obj.get_boundary_ops(boundary); [e_v,d_v,s_v,gamm_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + gamm_u = obj.getBoundaryBorrowing(boundary); + s_u = obj.getBoundarySign(boundary); + halfnorm_inv = obj.getHalfnormInv(boundary); + + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); + s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); + tuning = 1.1; alpha_u = obj.alpha; @@ -183,34 +195,109 @@ penalty = halfnorm_inv*(-tau*e_v' - sig*alpha_v*d_v'); end - % Ruturns the boundary ops and sign for the boundary specified by the string boundary. - % The right boundary is considered the positive boundary - function [e,d,s,gamm, halfnorm_inv] = get_boundary_ops(obj,boundary) + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string or a cell array of strings + % boundary -- string + function varargout = getBoundaryOperator(obj, op, boundary) + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'w' + e = obj.e_w; + case 'e' + e = obj.e_e; + case 's' + e = obj.e_s; + case 'n' + e = obj.e_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd' + switch boundary + case 'w' + d = obj.d1_w; + case 'e' + d = obj.d1_e; + case 's' + d = obj.d1_s; + case 'n' + d = obj.d1_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d; + end + end + + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) + switch boundary case 'w' - e = obj.e_w; - d = obj.d1_w; - s = -1; - gamm = obj.gamm_x; - halfnorm_inv = obj.Hix; + H_b = obj.H_y; case 'e' - e = obj.e_e; - d = obj.d1_e; - s = 1; - gamm = obj.gamm_x; - halfnorm_inv = obj.Hix; + H_b = obj.H_y; case 's' - e = obj.e_s; - d = obj.d1_s; - s = -1; - gamm = obj.gamm_y; - halfnorm_inv = obj.Hiy; + H_b = obj.H_x; case 'n' - e = obj.e_n; - d = obj.d1_n; + H_b = obj.H_x; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns borrowing constant gamma + % boundary -- string + function gamm = getBoundaryBorrowing(obj, boundary) + switch boundary + case {'w','e'} + gamm = obj.gamm_x; + case {'s','n'} + gamm = obj.gamm_y; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) + switch boundary + case {'e','n'} s = 1; - gamm = obj.gamm_y; - halfnorm_inv = obj.Hiy; + case {'w','s'} + s = -1; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns the halfnorm_inv used in SATs. TODO: better notation + function Hinv = getHalfnormInv(obj, boundary) + switch boundary + case 'w' + Hinv = obj.Hix; + case 'e' + Hinv = obj.Hix; + case 's' + Hinv = obj.Hiy; + case 'n' + Hinv = obj.Hiy; otherwise error('No such boundary: boundary = %s',boundary); end