Mercurial > repos > public > sbplib
changeset 1000:bd54cb25d96b feature/getBoundaryOp
Merge with default
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Mon, 14 Jan 2019 11:12:42 -0800 |
| parents | 337c4d1dcef5 (diff) a72038b1f709 (current diff) |
| children | 514a98f9f90d e512714fb890 |
| files | |
| diffstat | 16 files changed, 1265 insertions(+), 468 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+multiblock/DiffOp.m Mon Jan 14 11:12:42 2019 -0800 @@ -129,19 +129,20 @@ % Get a boundary operator specified by opName for the given boundary/BoundaryGroup function op = getBoundaryOperator(obj, opName, boundary) + blockmatrixDiv = obj.blockmatrixDiv{1}; + switch class(boundary) case 'cell' - localOpName = [opName '_' boundary{2}]; blockId = boundary{1}; - localOp = obj.diffOps{blockId}.(localOpName); + localOp = obj.diffOps{blockId}.getBoundaryOperator(opName, boundary{2}); - div = {obj.blockmatrixDiv{1}, size(localOp,2)}; + div = {blockmatrixDiv, size(localOp,2)}; blockOp = blockmatrix.zero(div); blockOp{blockId,1} = localOp; op = blockmatrix.toMatrix(blockOp); return case 'multiblock.BoundaryGroup' - op = sparse(size(obj.D,1),0); + op = sparse(sum(blockmatrixDiv),0); for i = 1:length(boundary) op = [op, obj.getBoundaryOperator(opName, boundary{i})]; end @@ -151,13 +152,10 @@ end function op = getBoundaryQuadrature(obj, boundary) - opName = 'H'; switch class(boundary) case 'cell' - localOpName = [opName '_' boundary{2}]; blockId = boundary{1}; - op = obj.diffOps{blockId}.(localOpName); - + op = obj.diffOps{blockId}.getBoundaryQuadrature(boundary{2}); return case 'multiblock.BoundaryGroup' N = length(boundary);
--- a/+scheme/Beam.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Beam.m Mon Jan 14 11:12:42 2019 -0800 @@ -86,7 +86,8 @@ function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dn'); - [e, d1, d2, d3, s] = obj.get_boundary_ops(boundary); + [e, d1, d2, d3] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary); + s = obj.getBoundarySign(boundary); gamm = obj.gamm; delt = obj.delt; @@ -173,14 +174,15 @@ 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] = obj.get_boundary_ops(boundary); - [e_v,d1_v,d2_v,d3_v,s_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); + [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary); + s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); alpha_u = obj.alpha; alpha_v = neighbour_scheme.alpha; - switch boundary case 'l' interface_opt = obj.opt.interface_l; @@ -234,22 +236,70 @@ penalty = -obj.Hi*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); 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, d1, d2, d3, s] = 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 ~ismember(boundary, {'l', 'r'}) + error('No such boundary: boundary = %s',boundary); + end + + if ~iscell(op) + op = {op}; + end + + for i = 1:numel(op) + switch op{i} + case 'e' + switch boundary + case 'l' + e = obj.e_l; + case 'r' + e = obj.e_r; + end + varargout{i} = e; + + case 'd1' + switch boundary + case 'l' + d1 = obj.d1_l; + case 'r' + d1 = obj.d1_r; + end + varargout{i} = d1; + end + + case 'd2' + switch boundary + case 'l' + d2 = obj.d2_l; + case 'r' + d2 = obj.d2_r; + end + varargout{i} = d2; + end + + case 'd3' + switch boundary + case 'l' + d3 = obj.d3_l; + case 'r' + d3 = obj.d3_r; + end + varargout{i} = d3; + end + end + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) switch boundary - case 'l' - e = obj.e_l; - d1 = obj.d1_l; - d2 = obj.d2_l; - d3 = obj.d3_l; + case {'r'} + s = 1; + case {'l'} s = -1; - case 'r' - e = obj.e_r; - d1 = obj.d1_r; - d2 = obj.d2_r; - d3 = obj.d3_r; - s = 1; otherwise error('No such boundary: boundary = %s',boundary); end
--- a/+scheme/Beam2d.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Beam2d.m Mon Jan 14 11:12:42 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/Elastic2dVariable.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Elastic2dVariable.m Mon Jan 14 11:12:42 2019 -0800 @@ -30,18 +30,10 @@ T_l, T_r tau_l, tau_r - H, Hi % Inner products - - phi % Borrowing constant for (d1 - e^T*D1) from R - gamma % Borrowing constant for d1 from M - H11 % First element of H + H, Hi, H_1D % Inner products + e_l, e_r - % Borrowing from H, M, and R - thH - thM - thR - e_l, e_r d1_l, d1_r % Normal derivatives at the boundary E % E{i}^T picks out component i @@ -50,22 +42,38 @@ % Kroneckered norms and coefficients RHOi_kron Hi_kron + + % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. + theta_R % Borrowing (d1- D1)^2 from R + theta_H % First entry in norm matrix + theta_M % Borrowing d1^2 from M. + + % Structures used for adjoint optimization + B end methods - function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) + % The coefficients can either be function handles or grid functions + function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet) default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); - default_arg('lambda_fun', @(x,y) 0*x+1); - default_arg('mu_fun', @(x,y) 0*x+1); - default_arg('rho_fun', @(x,y) 0*x+1); + default_arg('lambda', @(x,y) 0*x+1); + default_arg('mu', @(x,y) 0*x+1); + default_arg('rho', @(x,y) 0*x+1); dim = 2; assert(isa(g, 'grid.Cartesian')) - lambda = grid.evalOn(g, lambda_fun); - mu = grid.evalOn(g, mu_fun); - rho = grid.evalOn(g, rho_fun); + if isa(lambda, 'function_handle') + lambda = grid.evalOn(g, lambda); + end + if isa(mu, 'function_handle') + mu = grid.evalOn(g, mu); + end + if isa(rho, 'function_handle') + rho = grid.evalOn(g, rho); + end + m = g.size(); m_tot = g.N(); @@ -87,15 +95,9 @@ % Borrowing constants for i = 1:dim - beta = ops{i}.borrowing.R.delta_D; - obj.H11{i} = ops{i}.borrowing.H11; - obj.phi{i} = beta/obj.H11{i}; - obj.gamma{i} = ops{i}.borrowing.M.d1; - - % Better names - obj.thR{i} = ops{i}.borrowing.R.delta_D; - obj.thM{i} = ops{i}.borrowing.M.d1; - obj.thH{i} = ops{i}.borrowing.H11; + obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D; + obj.theta_H{i} = h(i)*ops{i}.borrowing.H11; + obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1; end I = cell(dim,1); @@ -183,6 +185,7 @@ obj.H_boundary = cell(dim,1); obj.H_boundary{1} = H{2}; obj.H_boundary{2} = H{1}; + obj.H_1D = {H{1}, H{2}}; % E{i}^T picks out component i. E = cell(dim,1); @@ -213,7 +216,7 @@ end end obj.D = D; - %=========================================% + %=========================================%' % Numerical traction operators for BC. % Because d1 =/= e0^T*D1, the numerical tractions are different @@ -237,20 +240,28 @@ tau_l{j} = cell(dim,1); tau_r{j} = cell(dim,1); + LAMBDA_l = e_l{j}'*LAMBDA*e_l{j}; + LAMBDA_r = e_r{j}'*LAMBDA*e_r{j}; + MU_l = e_l{j}'*MU*e_l{j}; + MU_r = e_r{j}'*MU*e_r{j}; + + [~, n_l] = size(e_l{j}); + [~, n_r] = size(e_r{j}); + % Loop over components for i = 1:dim - tau_l{j}{i} = sparse(m_tot,dim*m_tot); - tau_r{j}{i} = sparse(m_tot,dim*m_tot); + tau_l{j}{i} = sparse(n_l, dim*m_tot); + tau_r{j}{i} = sparse(n_r, dim*m_tot); for k = 1:dim T_l{j}{i,k} = ... - -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... - -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... - -d(i,k)*MU*e_l{j}*d1_l{j}'; + -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})... + -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})... + -d(i,k)*MU_l*d1_l{j}'; T_r{j}{i,k} = ... - d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... - +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... - +d(i,k)*MU*e_r{j}*d1_r{j}'; + d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})... + +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})... + +d(i,k)*MU_r*d1_r{j}'; tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; @@ -258,6 +269,19 @@ end end + + % Transpose T and tau to match boundary operator convention + for i = 1:dim + for j = 1:dim + tau_l{i}{j} = transpose(tau_l{i}{j}); + tau_r{i}{j} = transpose(tau_r{i}{j}); + for k = 1:dim + T_l{i}{j,k} = transpose(T_l{i}{j,k}); + T_r{i}{j,k} = transpose(T_r{i}{j,k}); + end + end + end + obj.T_l = T_l; obj.T_r = T_r; obj.tau_l = tau_l; @@ -275,6 +299,44 @@ obj.grid = g; obj.dim = dim; + % B, used for adjoint optimization + B = cell(dim, 1); + for i = 1:dim + B{i} = cell(m_tot, 1); + end + + for i = 1:dim + for j = 1:m_tot + B{i}{j} = sparse(m_tot, m_tot); + end + end + + ind = grid.funcToMatrix(g, 1:m_tot); + + % Direction 1 + for k = 1:m(1) + c = sparse(m(1),1); + c(k) = 1; + [~, B_1D] = ops{1}.D2(c); + for l = 1:m(2) + p = ind(:,l); + B{1}{(k-1)*m(2) + l}(p, p) = B_1D; + end + end + + % Direction 2 + for k = 1:m(2) + c = sparse(m(2),1); + c(k) = 1; + [~, B_1D] = ops{2}.D2(c); + for l = 1:m(1) + p = ind(l,:); + B{2}{(l-1)*m(2) + k}(p, p) = B_1D; + end + end + + obj.B = B; + end @@ -295,7 +357,8 @@ % j is the coordinate direction of the boundary j = obj.get_boundary_number(boundary); - [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); + [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); + E = obj.E; Hi = obj.Hi; @@ -316,33 +379,20 @@ % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - phi = obj.phi{j}; - h = obj.h(j); - h11 = obj.H11{j}*h; - gamma = obj.gamma{j}; - - a_lambda = dim/h11 + 1/(h11*phi); - a_mu_i = 2/(gamma*h); - a_mu_ij = 2/h11 + 1/(h11*phi); - - d = @kroneckerDelta; % Kronecker delta - db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta - alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... - + d(i,j)* a_mu_i*MU ... - + db(i,j)*a_mu_ij*MU ); + alpha = obj.getBoundaryOperator('alpha', boundary); % Loop over components that Dirichlet penalties end up on for i = 1:dim - C = T{k,i}; - A = -d(i,k)*alpha(i,j); - B = A + C; + C = transpose(T{k,i}); + A = -tuning*e*transpose(alpha{i,k}); + B = A + e*C; closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; end % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); + closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}'; penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; % Unknown boundary condition @@ -351,82 +401,59 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % type Struct that specifies the interface coupling. + % Fields: + % -- tuning: penalty strength, defaults to 1.2 + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.tuning = 1.2; + defaultType.interpolation = 'none'; + default_struct('type', defaultType); + + switch type.interpolation + case {'none', ''} + [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); + case {'op','OP'} + [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); + otherwise + error('Unknown type of interpolation: %s ', type.interpolation); + end + end + + function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) + tuning = type.tuning; + % u denotes the solution in the own domain % v denotes the solution in the neighbour domain % Operators without subscripts are from the own domain. - tuning = 1.2; - - % j is the coordinate direction of the boundary - j = obj.get_boundary_number(boundary); - j_v = neighbour_scheme.get_boundary_number(neighbour_boundary); % Get boundary operators - [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); - [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary); + [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary); + [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); % Operators and quantities that correspond to the own domain only - Hi = obj.Hi; - RHOi = obj.RHOi; - dim = obj.dim; - - %--- Other operators ---- - m_tot_u = obj.grid.N(); - E = obj.E; - LAMBDA_u = obj.LAMBDA; - MU_u = obj.MU; - lambda_u = e'*LAMBDA_u*e; - mu_u = e'*MU_u*e; + Hi = obj.Hi_kron; + RHOi = obj.RHOi_kron; - m_tot_v = neighbour_scheme.grid.N(); - E_v = neighbour_scheme.E; - LAMBDA_v = neighbour_scheme.LAMBDA; - MU_v = neighbour_scheme.MU; - lambda_v = e_v'*LAMBDA_v*e_v; - mu_v = e_v'*MU_v*e_v; - %------------------------- + % Penalty strength operators + alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary); + alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary); - % Borrowing constants - h_u = obj.h(j); - thR_u = obj.thR{j}*h_u; - thM_u = obj.thM{j}*h_u; - thH_u = obj.thH{j}*h_u; - - h_v = neighbour_scheme.h(j_v); - thR_v = neighbour_scheme.thR{j_v}*h_v; - thH_v = neighbour_scheme.thH{j_v}*h_v; - thM_v = neighbour_scheme.thM{j_v}*h_v; + closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e'); + penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v'); - % alpha = penalty strength for normal component, beta for tangential - alpha_u = dim*lambda_u/(4*thH_u) + lambda_u/(4*thR_u) + mu_u/(2*thM_u); - alpha_v = dim*lambda_v/(4*thH_v) + lambda_v/(4*thR_v) + mu_v/(2*thM_v); - beta_u = mu_u/(2*thH_u) + mu_u/(4*thR_u); - beta_v = mu_v/(2*thH_v) + mu_v/(4*thR_v); - alpha = alpha_u + alpha_v; - beta = beta_u + beta_v; - - d = @kroneckerDelta; % Kronecker delta - db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta - strength = @(i,j) tuning*(d(i,j)*alpha + db(i,j)*beta); + closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau'; + penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v'; - % Preallocate - closure = sparse(dim*m_tot_u, dim*m_tot_u); - penalty = sparse(dim*m_tot_u, dim*m_tot_v); - - % Loop over components that penalties end up on - for i = 1:dim - closure = closure - E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e'*E{i}'; - penalty = penalty + E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e_v'*E_v{i}'; + closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e'; + penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v'; - closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i}; - penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i}; + end - % Loop over components that we have interface conditions on - for k = 1:dim - closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}'; - penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}'; - end - end + function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) + error('Non-conforming interfaces not implemented yet.'); end % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. @@ -451,7 +478,8 @@ % Returns the boundary operator op for the boundary specified by the string boundary. % op: may be a cell array of strings - function [varargout] = get_boundary_operator(obj, op, boundary) + % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() + function [varargout] = getBoundaryOperator(obj, op, boundary) switch boundary case {'w','W','west','West', 'e', 'E', 'east', 'East'} @@ -466,45 +494,125 @@ op = {op}; end - for i = 1:length(op) - switch op{i} + for k = 1:length(op) + switch op{k} case 'e' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.e_l{j}; + varargout{k} = obj.e_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.e_r{j}; + varargout{k} = obj.e_r{j}; end + + case 'e_tot' + e = obj.getBoundaryOperator('e', boundary); + I_dim = speye(obj.dim, obj.dim); + varargout{k} = kron(e, I_dim); + case 'd' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.d1_l{j}; + varargout{k} = obj.d1_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.d1_r{j}; + varargout{k} = obj.d1_r{j}; end - case 'H' - varargout{i} = obj.H_boundary{j}; + case 'T' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.T_l{j}; + varargout{k} = obj.T_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.T_r{j}; + varargout{k} = obj.T_r{j}; end + case 'tau' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.tau_l{j}; + varargout{k} = obj.tau_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.tau_r{j}; + varargout{k} = obj.tau_r{j}; + end + + case 'tau_tot' + [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary); + + I_dim = speye(obj.dim, obj.dim); + e_tot = kron(e, I_dim); + E = obj.E; + tau_tot = (e_tot'*E{1}*e*tau{1}')'; + for i = 2:obj.dim + tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')'; end + varargout{k} = tau_tot; + + case 'H' + varargout{k} = obj.H_boundary{j}; + + case 'alpha' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + e = obj.getBoundaryOperator('e', boundary); + + LAMBDA = obj.LAMBDA; + MU = obj.MU; + + dim = obj.dim; + theta_R = obj.theta_R{j}; + theta_H = obj.theta_H{j}; + theta_M = obj.theta_M{j}; + + a_lambda = dim/theta_H + 1/theta_R; + a_mu_i = 2/theta_M; + a_mu_ij = 2/theta_H + 1/theta_R; + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + alpha = cell(obj.dim, obj.dim); + + alpha_func = @(i,j) d(i,j)* a_lambda*LAMBDA ... + + d(i,j)* a_mu_i*MU ... + + db(i,j)*a_mu_ij*MU; + for i = 1:obj.dim + for l = 1:obj.dim + alpha{i,l} = d(i,l)*alpha_func(i,j)*e; + end + end + + varargout{k} = alpha; + + case 'alpha_tot' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary); + E = obj.E; + [m, n] = size(alpha{1,1}); + alpha_tot = sparse(m*obj.dim, n*obj.dim); + for i = 1:obj.dim + for l = 1:obj.dim + alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; + end + end + varargout{k} = alpha_tot; + otherwise - error(['No such operator: operator = ' op{i}]); + error(['No such operator: operator = ' op{k}]); end end end + function H = 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; + otherwise + error('No such boundary: boundary = %s',boundary); + end + H = obj.H_boundary{j}; + I_dim = speye(obj.dim, obj.dim); + H = kron(H, I_dim); + end + function N = size(obj) N = obj.dim*prod(obj.m); end
--- a/+scheme/Euler1d.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Euler1d.m Mon Jan 14 11:12:42 2019 -0800 @@ -201,7 +201,8 @@ % Enforces the boundary conditions % w+ = R*w- + g(t) function closure = boundary_condition(obj,boundary, type, varargin) - [e_s,e_S,s] = obj.get_boundary_ops(boundary); + [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary); + s = obj.getBoundarySign(boundary); % Boundary condition on form % w_in = R*w_out + g, where g is data @@ -232,7 +233,8 @@ % % Returns closure(q,g) function closure = boundary_condition_L(obj, boundary, L_fun, p_in) - [e_s,e_S,s] = obj.get_boundary_ops(boundary); + [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary); + s = obj.getBoundarySign(boundary); p_ot = 1:3; p_ot(p_in) = []; @@ -273,7 +275,8 @@ % Return closure(q,g) function closure = boundary_condition_char(obj,boundary) - [e_s,e_S,s] = obj.get_boundary_ops(boundary); + [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary); + s = obj.getBoundarySign(boundary); function o = closure_fun(q, w_data) q_s = e_S'*q; @@ -314,7 +317,7 @@ % Return closure(q,[v; p]) function closure = boundary_condition_inflow(obj, boundary) - [~,~,s] = obj.get_boundary_ops(boundary); + s = obj.getBoundarySign(boundary); switch s case -1 @@ -335,7 +338,7 @@ % Return closure(q, p) function closure = boundary_condition_outflow(obj, boundary) - [~,~,s] = obj.get_boundary_ops(boundary); + s = obj.getBoundarySign(boundary); switch s case -1 @@ -352,7 +355,7 @@ % Return closure(q,[v; rho]) function closure = boundary_condition_inflow_rho(obj, boundary) - [~,~,s] = obj.get_boundary_ops(boundary); + s = obj.getBoundarySign(boundary); switch s case -1 @@ -372,7 +375,7 @@ % Return closure(q,rho) function closure = boundary_condition_outflow_rho(obj, boundary) - [~,~,s] = obj.get_boundary_ops(boundary); + s = obj.getBoundarySign(boundary); switch s case -1 @@ -388,7 +391,8 @@ % Set wall boundary condition v = 0. function closure = boundary_condition_wall(obj,boundary) - [e_s,e_S,s] = obj.get_boundary_ops(boundary); + [e_s, e_S] = obj.getBoundaryOperator({'e', 'E'}, boundary); + s = obj.getBoundarySign(boundary); % Vill vi sätta penalty på karateristikan som är nära noll också eller vill % vi låta den vara fri? @@ -478,18 +482,50 @@ 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,E,s] = 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 'l' + e = obj.e_l; + case 'r' + e = obj.e_r; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'E' + switch boundary + case 'l' + E = obj.e_L; + case 'r' + E = obj.e_R; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = E; + end + end + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) switch boundary - case 'l' - e = obj.e_l; - E = obj.e_L; + case {'r'} + s = 1; + case {'l'} s = -1; - case 'r' - e = obj.e_r; - E = obj.e_R; - s = 1; otherwise error('No such boundary: boundary = %s',boundary); end
--- a/+scheme/Heat2dCurvilinear.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Heat2dCurvilinear.m Mon Jan 14 11:12:42 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 15:00:12 2019 +0100 +++ b/+scheme/Heat2dVariable.m Mon Jan 14 11:12:42 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 15:00:12 2019 +0100 +++ b/+scheme/Hypsyst2d.m Mon Jan 14 11:12:42 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 15:00:12 2019 +0100 +++ b/+scheme/Hypsyst2dCurve.m Mon Jan 14 11:12:42 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/Laplace1d.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Laplace1d.m Mon Jan 14 11:12:42 2019 -0800 @@ -56,7 +56,8 @@ default_arg('type','neumann'); default_arg('data',0); - [e,d,s] = obj.get_boundary_ops(boundary); + [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + s = obj.getBoundarySign(boundary); switch type % Dirichlet boundary condition @@ -86,10 +87,11 @@ 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, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + s_u = obj.getBoundarySign(boundary); - [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); - [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); a_u = obj.a; a_v = neighbour_scheme.a; @@ -111,18 +113,50 @@ penalty = obj.Hi*(-tau*e_v' + sig*a_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] = 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 'l' + e = obj.e_l; + case 'r' + e = obj.e_r; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd' + switch boundary + case 'l' + d = obj.d_l; + case 'r' + d = obj.d_r; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d; + end + end + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) switch boundary - case 'l' - e = obj.e_l; - d = obj.d_l; + case {'r'} + s = 1; + case {'l'} s = -1; - case 'r' - e = obj.e_r; - d = obj.d_r; - s = 1; otherwise error('No such boundary: boundary = %s',boundary); end
--- a/+scheme/LaplaceCurvilinear.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/LaplaceCurvilinear.m Mon Jan 14 11:12:42 2019 -0800 @@ -238,7 +238,9 @@ default_arg('type','neumann'); default_arg('parameter', []); - [e, d, gamm, H_b, ~] = obj.get_boundary_ops(boundary); + [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_b = obj.getBoundaryQuadrature(boundary); + gamm = obj.getBoundaryBorrowing(boundary); switch type % Dirichlet boundary condition case {'D','d','dirichlet'} @@ -298,8 +300,15 @@ % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); - [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_b_u = obj.getBoundaryQuadrature(boundary); + I_u = obj.getBoundaryIndices(boundary); + gamm_u = obj.getBoundaryBorrowing(boundary); + + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); + I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); + gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); u = obj; v = neighbour_scheme; @@ -336,8 +345,16 @@ % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - [e_u, d_u, gamm_u, H_b_u, I_u] = obj.get_boundary_ops(boundary); - [e_v, d_v, gamm_v, H_b_v, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + H_b_u = obj.getBoundaryQuadrature(boundary); + I_u = obj.getBoundaryIndices(boundary); + gamm_u = obj.getBoundaryBorrowing(boundary); + + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); + I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); + gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary); + % Find the number of grid points along the interface m_u = size(e_u, 2); @@ -378,43 +395,98 @@ end - % Returns the boundary ops and sign for the boundary specified by the string boundary. - % The right boundary is considered the positive 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 % - % I -- the indices of the boundary points in the grid matrix - function [e, d, gamm, H_b, I] = get_boundary_ops(obj, boundary) - ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); + % boundary -- string + function H_b = getBoundaryQuadrature(obj, boundary) switch boundary case 'w' - e = obj.e_w; - d = obj.d_w; H_b = obj.H_w; + case 'e' + H_b = obj.H_e; + case 's' + H_b = obj.H_s; + case 'n' + H_b = obj.H_n; + otherwise + error('No such boundary: boundary = %s',boundary); + end + end + + % Returns the indices of the boundary points in the grid matrix + % boundary -- string + function I = getBoundaryIndices(obj, boundary) + ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); + switch boundary + case 'w' I = ind(1,:); case 'e' - e = obj.e_e; - d = obj.d_e; - H_b = obj.H_e; I = ind(end,:); case 's' - e = obj.e_s; - d = obj.d_s; - H_b = obj.H_s; I = ind(:,1)'; case 'n' - e = obj.e_n; - d = obj.d_n; - H_b = obj.H_n; I = ind(:,end)'; 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_u; case {'s','n'} gamm = obj.gamm_v; + otherwise + error('No such boundary: boundary = %s',boundary); end end
--- a/+scheme/Schrodinger.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Schrodinger.m Mon Jan 14 11:12:42 2019 -0800 @@ -67,7 +67,8 @@ default_arg('type','dirichlet'); default_arg('data',0); - [e,d,s] = obj.get_boundary_ops(boundary); + [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + s = obj.getBoundarySign(boundary); switch type % Dirichlet boundary condition @@ -93,8 +94,11 @@ 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,d_u,s_u] = obj.get_boundary_ops(boundary); - [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); + [e_u, d_u] = obj.getBoundaryOperator({'e', 'd'}, boundary); + s_u = obj.getBoundarySign(boundary); + + [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); a = -s_u* 1/2 * 1i ; b = a'; @@ -106,18 +110,50 @@ penalty = obj.Hi * (-tau*e_v' - sig*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] = 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 'l' + e = obj.e_l; + case 'r' + e = obj.e_r; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + + case 'd' + switch boundary + case 'l' + d = obj.d1_l; + case 'r' + d = obj.d1_r; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = d; + end + end + end + + % Returns the boundary sign. The right boundary is considered the positive boundary + % boundary -- string + function s = getBoundarySign(obj, boundary) switch boundary - case 'l' - e = obj.e_l; - d = obj.d1_l; + case {'r'} + s = 1; + case {'l'} s = -1; - case 'r' - e = obj.e_r; - d = obj.d1_r; - s = 1; otherwise error('No such boundary: boundary = %s',boundary); end
--- a/+scheme/Schrodinger2d.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Schrodinger2d.m Mon Jan 14 11:12:42 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/Utux.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Utux.m Mon Jan 14 11:12:42 2019 -0800 @@ -72,6 +72,31 @@ 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 'l' + e = obj.e_l; + case 'r' + e = obj.e_r; + otherwise + error('No such boundary: boundary = %s',boundary); + end + varargout{i} = e; + end + end + end + function N = size(obj) N = obj.m; end
--- a/+scheme/Utux2d.m Tue Jan 08 15:00:12 2019 +0100 +++ b/+scheme/Utux2d.m Mon Jan 14 11:12:42 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 15:00:12 2019 +0100 +++ b/+scheme/Wave2d.m Mon Jan 14 11:12:42 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