Mercurial > repos > public > sbplib
changeset 1108:5ec23b9bf360 feature/laplace_curvilinear_test
Merge with default
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 10 Apr 2019 11:00:27 -0700 |
| parents | 867307f4d80f (current diff) 00203fcc962f (diff) |
| children | 01d28cfafe7c |
| files | +scheme/Beam2d.m +scheme/LaplaceCurvilinear.m +scheme/TODO.txt +scheme/Wave2dCurve.m +scheme/error1d.m +scheme/error2d.m +scheme/errorMax.m +scheme/errorRelative.m +scheme/errorSbp.m +scheme/errorVector.m |
| diffstat | 31 files changed, 608 insertions(+), 1170 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/+domain/Line.m Wed Apr 10 11:00:27 2019 -0700 @@ -0,0 +1,153 @@ +classdef Line < multiblock.Definition + properties + + xlims + blockNames % Cell array of block labels + nBlocks + connections % Cell array specifying connections between blocks + boundaryGroups % Structure of boundaryGroups + + end + + + methods + % Creates a divided line + % x is a vector of boundary and interface positions. + % blockNames: cell array of labels. The id is default. + function obj = Line(x,blockNames) + default_arg('blockNames',[]); + + N = length(x)-1; % number of blocks in the x direction. + + if ~issorted(x) + error('The elements of x seem to be in the wrong order'); + end + + % Dimensions of blocks and number of points + blockTi = cell(N,1); + xlims = cell(N,1); + for i = 1:N + xlims{i} = {x(i), x(i+1)}; + end + + % Interface couplings + conn = cell(N,N); + for i = 1:N + conn{i,i+1} = {'r','l'}; + end + + % Block names (id number as default) + if isempty(blockNames) + obj.blockNames = cell(1, N); + for i = 1:N + obj.blockNames{i} = sprintf('%d', i); + end + else + assert(length(blockNames) == N); + obj.blockNames = blockNames; + end + nBlocks = N; + + % Boundary groups + boundaryGroups = struct(); + L = { {1, 'l'} }; + R = { {N, 'r'} }; + boundaryGroups.L = multiblock.BoundaryGroup(L); + boundaryGroups.R = multiblock.BoundaryGroup(R); + boundaryGroups.all = multiblock.BoundaryGroup([L,R]); + + obj.connections = conn; + obj.nBlocks = nBlocks; + obj.boundaryGroups = boundaryGroups; + obj.xlims = xlims; + + end + + + % Returns a multiblock.Grid given some parameters + % ms: cell array of m values + % For same m in every block, just input one scalar. + function g = getGrid(obj, ms, varargin) + + default_arg('ms',21) + + % Extend ms if input is a single scalar + if (numel(ms) == 1) && ~iscell(ms) + m = ms; + ms = cell(1,obj.nBlocks); + for i = 1:obj.nBlocks + ms{i} = m; + end + end + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + grids{i} = grid.equidistant(ms{i}, obj.xlims{i}); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + % Returns a multiblock.Grid given some parameters + % ms: cell array of m values + % For same m in every block, just input one scalar. + function g = getStaggeredGrid(obj, ms, varargin) + + default_arg('ms',21) + + % Extend ms if input is a single scalar + if (numel(ms) == 1) && ~iscell(ms) + m = ms; + ms = cell(1,obj.nBlocks); + for i = 1:obj.nBlocks + ms{i} = m; + end + end + + grids = cell(1, obj.nBlocks); + for i = 1:obj.nBlocks + [g_primal, g_dual] = grid.primalDual1D(ms{i}, obj.xlims{i}); + grids{i} = grid.Staggered1d(g_primal, g_dual); + end + + g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); + end + + % label is the type of label used for plotting, + % default is block name, 'id' show the index for each block. + function show(obj, label) + default_arg('label', 'name') + + m = 10; + figure + for i = 1:obj.nBlocks + x = linspace(obj.xlims{i}{1}, obj.xlims{i}{2}, m); + y = 0*x + 0.05* ( (-1)^i + 1 ) ; + plot(x,y,'+'); + hold on + end + hold off + + switch label + case 'name' + labels = obj.blockNames; + case 'id' + labels = {}; + for i = 1:obj.nBlocks + labels{i} = num2str(i); + end + otherwise + axis equal + return + end + + legend(labels) + axis equal + end + + % Returns the grid size of each block in a cell array + % The input parameters are determined by the subclass + function ms = getGridSizes(obj, varargin) + end + end +end
--- a/+multiblock/DiffOp.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+multiblock/DiffOp.m Wed Apr 10 11:00:27 2019 -0700 @@ -129,20 +129,19 @@ % 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' blockId = boundary{1}; localOp = obj.diffOps{blockId}.getBoundaryOperator(opName, boundary{2}); - div = {blockmatrixDiv, size(localOp,2)}; + div = {obj.blockmatrixDiv{1}, size(localOp,2)}; blockOp = blockmatrix.zero(div); blockOp{blockId,1} = localOp; op = blockmatrix.toMatrix(blockOp); return case 'multiblock.BoundaryGroup' - op = sparse(sum(blockmatrixDiv),0); + op = sparse(size(obj.D,1),0); for i = 1:length(boundary) op = [op, obj.getBoundaryOperator(opName, boundary{i})]; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/LaplaceSquared.m Wed Apr 10 11:00:27 2019 -0700 @@ -0,0 +1,105 @@ +classdef LaplaceSquared < scheme.Scheme + properties + grid + order + laplaceDiffOp + + D + H + Hi + + a,b + end + + methods + % Discretisation of a*nabla*b*nabla + function obj = LaplaceSquared(g, order, a, b, opGen) + default_arg('order', 4); + default_arg('a', 1); + default_arg('b', 1); + default_arg('opGen', @sbp.D4Variable); + + if isscalar(a) + a = grid.evalOn(g, a); + end + + if isscalar(b) + b = grid.evalOn(g, b); + end + + obj.grid = g; + obj.order = order; + obj.a = a; + obj.b = b; + + obj.laplaceDiffOp = multiblock.Laplace(g, order, 1, 1, opGen); + + obj.H = obj.laplaceDiffOp.H; + obj.Hi = spdiag(1./diag(obj.H)); + + A = spdiag(a); + B = spdiag(b); + + D_laplace = obj.laplaceDiffOp.D; + obj.D = A*D_laplace*B*D_laplace; + end + + function s = size(obj) + s = size(obj.laplaceDiffOp); + end + + function op = getBoundaryOperator(obj, opName, boundary) + switch opName + case 'e' + op = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary); + case 'd1' + op = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary); + case 'd2' + e = getBoundaryOperator(obj.laplaceDiffOp, 'e', boundary); + op = (e'*obj.laplaceDiffOp.D)'; + case 'd3' + d1 = getBoundaryOperator(obj.laplaceDiffOp, 'd', boundary); + op = (d1'*spdiag(obj.b)*obj.laplaceDiffOp.D)'; + end + end + + function op = getBoundaryQuadrature(obj, boundary) + op = getBoundaryQuadrature(obj.laplaceDiffOp, boundary); + end + + function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition + switch type + case 'e' + error('Bc of type ''e'' not implemented') + case 'd1' + error('Bc of type ''d1'' not implemented') + case 'd2' + e = obj.getBoundaryOperator('e', boundary); + d1 = obj.getBoundaryOperator('d1', boundary); + d2 = obj.getBoundaryOperator('d2', boundary); + H_b = obj.getBoundaryQuadrature(boundary); + + A = spdiag(obj.a); + B_b = spdiag(e'*obj.b); + + tau = obj.Hi*A*d1*B_b*H_b; + closure = tau*d2'; + penalty = -tau; + case 'd3' + e = obj.getBoundaryOperator('e', boundary); + d3 = obj.getBoundaryOperator('d3', boundary); + H_b = obj.getBoundaryQuadrature(boundary); + + A = spdiag(obj.a); + + tau = -obj.Hi*A*e*H_b; + closure = tau*d3'; + penalty = -tau; + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + error('Not implemented') + end + end +end
--- a/+parametrization/Curve.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+parametrization/Curve.m Wed Apr 10 11:00:27 2019 -0700 @@ -103,7 +103,10 @@ % Construct arcLength function using splines tvec = linspace(0,1,N); arcVec = obj.arcLength(0,tvec); - tFunc = spline(arcVec,tvec); % t as a function of arcLength + + % t as a function of arcLength. Monotonicity-preserving cubic splines. + tFunc = @(arcLen) pchip(arcVec,tvec,arcLen); + L = obj.arcLength(0,1); arcPar = @(s) tFunc(s*L); @@ -349,8 +352,6 @@ end end - - function g_norm = normalize(g0) g1 = g0(1,:); g2 = g0(2,:);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+parametrization/dataSpline.m Wed Apr 10 11:00:27 2019 -0700 @@ -0,0 +1,19 @@ +% dataSpline calculates a Curve through the points f_i using cubic spline interpolation. +% The spline curve is parametrized with the arc length parametrization +% to facilitate better grids. +% +% f - m x D matrix of m points in D dimensions +function C = dataSpline(f) + m = size(f, 1); + + t = linspace(0,1,m); + + pp_g = spapi(4, t, f'); + pp_gp = fnder(pp_g); + + g = @(t) fnval(pp_g, t); + gp = @(t) fnval(pp_gp, t); + + C = parametrization.Curve(g, gp); + C = C.arcLengthParametrization(); +end
--- a/+scheme/Beam.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Beam.m Wed Apr 10 11:00:27 2019 -0700 @@ -86,7 +86,10 @@ function [closure, penalty] = boundary_condition(obj,boundary,type) default_arg('type','dn'); - [e, d1, d2, d3] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary); + e = obj.getBoundaryOperator('e', boundary); + d1 = obj.getBoundaryOperator('d1', boundary); + d2 = obj.getBoundaryOperator('d2', boundary); + d3 = obj.getBoundaryOperator('d3', boundary); s = obj.getBoundarySign(boundary); gamm = obj.gamm; delt = obj.delt; @@ -125,7 +128,7 @@ closure = obj.Hi*(tau*d2' + sig*d3'); penalty{1} = -obj.Hi*tau; - penalty{1} = -obj.Hi*sig; + penalty{2} = -obj.Hi*sig; case 'e' alpha = obj.alpha; @@ -174,10 +177,16 @@ 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] = obj.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, boundary); + e_u = obj.getBoundaryOperator('e', boundary); + d1_u = obj.getBoundaryOperator('d1', boundary); + d2_u = obj.getBoundaryOperator('d2', boundary); + d3_u = obj.getBoundaryOperator('d3', boundary); s_u = obj.getBoundarySign(boundary); - [e_v, d1_v, d2_v, d3_v] = neighbour_scheme.getBoundaryOperator({'e', 'd1', 'd2', 'd3'}, neighbour_boundary); + e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); + d1_v = neighbour_scheme.getBoundaryOperator('d1', neighbour_boundary); + d2_v = neighbour_scheme.getBoundaryOperator('d2', neighbour_boundary); + d3_v = neighbour_scheme.getBoundaryOperator('d3', neighbour_boundary); s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); alpha_u = obj.alpha; @@ -237,71 +246,36 @@ end % Returns the boundary operator op for the boundary specified by the string boundary. - % op -- string or a cell array of strings + % op -- string % boundary -- string - function varargout = getBoundaryOperator(obj, op, boundary) - - if ~ismember(boundary, {'l', 'r'}) - error('No such boundary: boundary = %s',boundary); - end + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e', 'd1', 'd2', 'd3'}) + assertIsMember(boundary, {'l', 'r'}) - 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; + o = obj.([op, '_', boundary]); + end - case 'd1' - switch boundary - case 'l' - d1 = obj.d1_l; - case 'r' - d1 = obj.d1_r; - end - varargout{i} = d1; - end + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + % Note: for 1d diffOps, the boundary quadrature is the scalar 1. + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) - 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 + H_b = 1; end % Returns the boundary sign. The right boundary is considered the positive boundary % boundary -- string function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + switch boundary case {'r'} s = 1; case {'l'} s = -1; - otherwise - error('No such boundary: boundary = %s',boundary); end end
--- a/+scheme/Beam2d.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,349 +0,0 @@ -classdef Beam2d < scheme.Scheme - properties - grid - order % Order accuracy for the approximation - - D % non-stabalized scheme operator - M % Derivative norm - alpha - - H % Discrete norm - Hi - H_x, H_y % Norms in the x and y directions - Hx,Hy % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. - Hi_x, Hi_y - Hix, Hiy - e_w, e_e, e_s, e_n - d1_w, d1_e, d1_s, d1_n - d2_w, d2_e, d2_s, d2_n - d3_w, d3_e, d3_s, d3_n - gamm_x, gamm_y - delt_x, delt_y - end - - methods - function obj = Beam2d(m,lim,order,alpha,opsGen) - default_arg('alpha',1); - default_arg('opsGen',@sbp.Higher); - - if ~isa(grid, 'grid.Cartesian') || grid.D() ~= 2 - error('Grid must be 2d cartesian'); - end - - obj.grid = grid; - obj.alpha = alpha; - obj.order = order; - - m_x = grid.m(1); - m_y = grid.m(2); - - h = grid.scaling(); - h_x = h(1); - h_y = h(2); - - ops_x = opsGen(m_x,h_x,order); - ops_y = opsGen(m_y,h_y,order); - - I_x = speye(m_x); - I_y = speye(m_y); - - D4_x = sparse(ops_x.derivatives.D4); - H_x = sparse(ops_x.norms.H); - Hi_x = sparse(ops_x.norms.HI); - e_l_x = sparse(ops_x.boundary.e_1); - e_r_x = sparse(ops_x.boundary.e_m); - d1_l_x = sparse(ops_x.boundary.S_1); - d1_r_x = sparse(ops_x.boundary.S_m); - d2_l_x = sparse(ops_x.boundary.S2_1); - d2_r_x = sparse(ops_x.boundary.S2_m); - d3_l_x = sparse(ops_x.boundary.S3_1); - d3_r_x = sparse(ops_x.boundary.S3_m); - - D4_y = sparse(ops_y.derivatives.D4); - H_y = sparse(ops_y.norms.H); - Hi_y = sparse(ops_y.norms.HI); - e_l_y = sparse(ops_y.boundary.e_1); - e_r_y = sparse(ops_y.boundary.e_m); - d1_l_y = sparse(ops_y.boundary.S_1); - d1_r_y = sparse(ops_y.boundary.S_m); - d2_l_y = sparse(ops_y.boundary.S2_1); - d2_r_y = sparse(ops_y.boundary.S2_m); - d3_l_y = sparse(ops_y.boundary.S3_1); - d3_r_y = sparse(ops_y.boundary.S3_m); - - - D4 = kr(D4_x, I_y) + kr(I_x, D4_y); - - % Norms - obj.H = kr(H_x,H_y); - obj.Hx = kr(H_x,I_x); - obj.Hy = kr(I_x,H_y); - obj.Hix = kr(Hi_x,I_y); - obj.Hiy = kr(I_x,Hi_y); - obj.Hi = kr(Hi_x,Hi_y); - - % Boundary operators - obj.e_w = kr(e_l_x,I_y); - obj.e_e = kr(e_r_x,I_y); - obj.e_s = kr(I_x,e_l_y); - obj.e_n = kr(I_x,e_r_y); - obj.d1_w = kr(d1_l_x,I_y); - obj.d1_e = kr(d1_r_x,I_y); - obj.d1_s = kr(I_x,d1_l_y); - obj.d1_n = kr(I_x,d1_r_y); - obj.d2_w = kr(d2_l_x,I_y); - obj.d2_e = kr(d2_r_x,I_y); - obj.d2_s = kr(I_x,d2_l_y); - obj.d2_n = kr(I_x,d2_r_y); - obj.d3_w = kr(d3_l_x,I_y); - obj.d3_e = kr(d3_r_x,I_y); - obj.d3_s = kr(I_x,d3_l_y); - obj.d3_n = kr(I_x,d3_r_y); - - obj.D = alpha*D4; - - obj.gamm_x = h_x*ops_x.borrowing.N.S2/2; - obj.delt_x = h_x^3*ops_x.borrowing.N.S3/2; - - obj.gamm_y = h_y*ops_y.borrowing.N.S2/2; - obj.delt_y = h_y^3*ops_y.borrowing.N.S3/2; - end - - - % Closure functions return the opertors applied to the own doamin to close the boundary - % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. - % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. - % data is a function returning the data that should be applied at the boundary. - % neighbour_scheme is an instance of Scheme that should be interfaced to. - % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty_e,penalty_d] = boundary_condition(obj,boundary,type,data) - default_arg('type','dn'); - default_arg('data',0); - - [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 - case {'dn'} - alpha = obj.alpha; - - % tau1 < -alpha^2/gamma - tuning = 1.1; - - tau1 = tuning * alpha/delt; - tau4 = s*alpha; - - sig2 = tuning * alpha/gamm; - sig3 = -s*alpha; - - tau = tau1*e+tau4*d3; - sig = sig2*d1+sig3*d2; - - closure = halfnorm_inv*(tau*e' + sig*d1'); - - pp_e = halfnorm_inv*tau; - pp_d = halfnorm_inv*sig; - switch class(data) - case 'double' - penalty_e = pp_e*data; - penalty_d = pp_d*data; - case 'function_handle' - penalty_e = @(t)pp_e*data(t); - penalty_d = @(t)pp_d*data(t); - otherwise - error('Wierd data argument!') - end - - % Unknown, boundary condition - otherwise - error('No such boundary condition: type = %s',type); - end - end - - 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] = 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; - - alpha_u = obj.alpha; - alpha_v = neighbour_scheme.alpha; - - tau1 = ((alpha_u/2)/delt_u + (alpha_v/2)/delt_v)/2*tuning; - % tau1 = (alpha_u/2 + alpha_v/2)/(2*delt_u)*tuning; - tau4 = s_u*alpha_u/2; - - sig2 = ((alpha_u/2)/gamm_u + (alpha_v/2)/gamm_v)/2*tuning; - sig3 = -s_u*alpha_u/2; - - phi2 = s_u*1/2; - - psi1 = -s_u*1/2; - - tau = tau1*e_u + tau4*d3_u; - sig = sig2*d1_u + sig3*d2_u ; - phi = phi2*d1_u ; - psi = psi1*e_u ; - - closure = halfnorm_inv*(tau*e_u' + sig*d1_u' + phi*alpha_u*d2_u' + psi*alpha_u*d3_u'); - penalty = -halfnorm_inv*(tau*e_v' + sig*d1_v' + phi*alpha_v*d2_v' + psi*alpha_v*d3_v'); - 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; - - 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' - 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; - case {'s','n'} - gamm = obj.gamm_y; - delt = obj.delt_y; - otherwise - error('No such boundary: boundary = %s',boundary); - end - end - - function N = size(obj) - N = prod(obj.m); - end - - end -end
--- a/+scheme/Elastic2dCurvilinear.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Elastic2dCurvilinear.m Wed Apr 10 11:00:27 2019 -0700 @@ -3,12 +3,12 @@ % Discretizes the elastic wave equation in curvilinear coordinates. % % Untransformed equation: -% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i +% rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i % % Transformed equation: -% J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j -% + dk J b_jk mu b_il dl u_j -% + dk J b_jk mu b_jl dl u_i +% J*rho u_{i,tt} = dk J b_ik lambda b_jl dl u_j +% + dk J b_jk mu b_il dl u_j +% + dk J b_jk mu b_jl dl u_i % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -49,7 +49,7 @@ e_l, e_r d1_l, d1_r % Normal derivatives at the boundary E % E{i}^T picks out component i - + H_boundary_l, H_boundary_r % Boundary inner products % Kroneckered norms and coefficients @@ -145,7 +145,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; @@ -327,12 +327,12 @@ for m = 1:dim for l = 1:dim - T_l{j}{i,k} = T_l{j}{i,k} + ... + T_l{j}{i,k} = T_l{j}{i,k} + ... -d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ... -d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}) ... -d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_l{m}*d1_l{m}' + db(m,j)*D1{m}); - T_r{j}{i,k} = T_r{j}{i,k} + ... + T_r{j}{i,k} = T_r{j}{i,k} + ... d(k,l)* J*b{i,j}*b{k,m}*LAMBDA*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ... d(k,l)* J*b{k,j}*b{i,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}) + ... d(i,k)* J*b{l,j}*b{l,m}*MU*(d(m,j)*e_r{m}*d1_r{m}' + db(m,j)*D1{m}); @@ -340,7 +340,7 @@ end T_l{j}{i,k} = inv(beta{j})*T_l{j}{i,k}; - T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k}; + T_r{j}{i,k} = inv(beta{j})*T_r{j}{i,k}; 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}'; @@ -387,7 +387,7 @@ % 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; @@ -423,20 +423,20 @@ 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 ); + + db(i,j)*a_mu_ij*MU ); % 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; - closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' ); + closure = closure + E{i}*RHOi*Hi*Ji*B'*e*H_gamma*(e'*E{k}' ); penalty = penalty - E{i}*RHOi*Hi*Ji*B'*e*H_gamma; - end + end % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} ); + closure = closure - E{k}*RHOi*Ji*Hi*e*H_gamma* (e'*tau{k} ); penalty = penalty + E{k}*RHOi*Ji*Hi*e*H_gamma; % Unknown boundary condition @@ -457,14 +457,14 @@ 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, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); + [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_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; @@ -480,7 +480,7 @@ lambda_v = e_v'*LAMBDA_v*e_v; mu_v = e_v'*MU_v*e_v; %------------------------- - + % Borrowing constants phi_u = obj.phi{j}; h_u = obj.h(j); @@ -493,7 +493,7 @@ gamma_v = neighbour_scheme.gamma{j_v}; % E > sum_i 1/(2*alpha_ij)*(tau_i)^2 - function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) + function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) th1 = h11/(2*dim); th2 = h11*phi/2; th3 = h*gamma; @@ -505,7 +505,7 @@ end [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u); - [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v); + [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v); sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4; sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4; @@ -527,9 +527,9 @@ % 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 + 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 end @@ -555,7 +555,7 @@ % 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) + function [varargout] = getBoundaryOperator(obj, op, boundary) switch boundary case {'w','W','west','West', 'e', 'E', 'east', 'East'} @@ -587,7 +587,7 @@ varargout{i} = obj.d1_r{j}; end case 'H' - switch boundary + switch boundary case {'w','W','west','West','s','S','south','South'} varargout{i} = obj.H_boundary_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} @@ -606,7 +606,7 @@ varargout{i} = obj.tau_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} varargout{i} = obj.tau_r{j}; - end + end otherwise error(['No such operator: operator = ' op{i}]); end @@ -614,6 +614,27 @@ end + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + function H = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + switch boundary + case {'w'} + H = H_boundary_l{1}; + case 'e' + H = H_boundary_r{1}; + case 's' + H = H_boundary_l{2}; + case 'n' + H = H_boundary_r{2}; + end + 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/Elastic2dVariable.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Elastic2dVariable.m Wed Apr 10 11:00:27 2019 -0700 @@ -429,8 +429,12 @@ % Operators without subscripts are from the own domain. % Get boundary operators - [e, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary); - [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary); + e = obj.getBoundaryOperator('e_tot', boundary); + tau = obj.getBoundaryOperator('tau_tot', boundary); + + e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary); + tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary); + H_gamma = obj.getBoundaryQuadrature(boundary); % Operators and quantities that correspond to the own domain only @@ -458,155 +462,149 @@ % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. function [j, nj] = get_boundary_number(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} + case {'w', 'e'} j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} + case {'s', 'n'} j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); end switch boundary - case {'w','W','west','West','s','S','south','South'} + case {'w', 's'} nj = -1; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + case {'e', 'n'} nj = 1; end end % Returns the boundary operator op for the boundary specified by the string boundary. - % op: may be a cell array of strings + % op -- string % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() function [varargout] = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'}) switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} + case {'w', 'e'} j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} + case {'s', 'n'} j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); - end - - if ~iscell(op) - op = {op}; end - for k = 1:length(op) - switch op{k} - case 'e' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.e_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.e_r{j}; - end + switch op + case 'e' + switch boundary + case {'w', 's'} + o = obj.e_l{j}; + case {'e', 'n'} + o = 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{k} = obj.d1_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.d1_r{j}; - end + case 'e_tot' + e = obj.getBoundaryOperator('e', boundary); + I_dim = speye(obj.dim, obj.dim); + o = kron(e, I_dim); - case 'T' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.T_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.T_r{j}; - end + case 'd' + switch boundary + case {'w', 's'} + o = obj.d1_l{j}; + case {'e', 'n'} + o = obj.d1_r{j}; + end - case 'tau' - switch boundary - case {'w','W','west','West','s','S','south','South'} - varargout{k} = obj.tau_l{j}; - case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{k} = obj.tau_r{j}; - end + case 'T' + switch boundary + case {'w', 's'} + o = obj.T_l{j}; + case {'e', 'n'} + o = obj.T_r{j}; + end - case 'tau_tot' - [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary); + case 'tau' + switch boundary + case {'w', 's'} + o = obj.tau_l{j}; + case {'e', 'n'} + o = obj.tau_r{j}; + end - 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 'tau_tot' + [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary); - case 'H' - varargout{k} = obj.H_boundary{j}; + 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 + o = tau_tot; - case 'alpha' - % alpha = alpha(i,j) is the penalty strength for displacement BC. - e = obj.getBoundaryOperator('e', boundary); - - LAMBDA = obj.LAMBDA; - MU = obj.MU; + case 'H' + o = obj.H_boundary{j}; - dim = obj.dim; - theta_R = obj.theta_R{j}; - theta_H = obj.theta_H{j}; - theta_M = obj.theta_M{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; - a_lambda = dim/theta_H + 1/theta_R; - a_mu_i = 2/theta_M; - a_mu_ij = 2/theta_H + 1/theta_R; + dim = obj.dim; + theta_R = obj.theta_R{j}; + theta_H = obj.theta_H{j}; + theta_M = obj.theta_M{j}; - d = @kroneckerDelta; % Kronecker delta - db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta - alpha = cell(obj.dim, obj.dim); + a_lambda = dim/theta_H + 1/theta_R; + a_mu_i = 2/theta_M; + a_mu_ij = 2/theta_H + 1/theta_R; - 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 + 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 - - varargout{k} = alpha; + end - 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 + o = 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 - varargout{k} = alpha_tot; - - otherwise - error(['No such operator: operator = ' op{k}]); - end + end + o = alpha_tot; end end + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string function H = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary - case {'w','W','west','West', 'e', 'E', 'east', 'East'} + case {'w','e'} j = 1; - case {'s','S','south','South', 'n', 'N', 'north', 'North'} + case {'s','n'} j = 2; - otherwise - error('No such boundary: boundary = %s',boundary); end H = obj.H_boundary{j}; I_dim = speye(obj.dim, obj.dim);
--- a/+scheme/Euler1d.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Euler1d.m Wed Apr 10 11:00:27 2019 -0700 @@ -518,6 +518,17 @@ end end + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + % Note: for 1d diffOps, the boundary quadrature is the scalar 1. + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + + H_b = 1; + end + % Returns the boundary sign. The right boundary is considered the positive boundary % boundary -- string function s = getBoundarySign(obj, boundary)
--- a/+scheme/Heat2dCurvilinear.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Heat2dCurvilinear.m Wed Apr 10 11:00:27 2019 -0700 @@ -321,6 +321,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) if ~iscell(op) op = {op}; @@ -338,8 +339,6 @@ e = obj.e_l{2}; case 'n' e = obj.e_r{2}; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; @@ -353,8 +352,6 @@ d = obj.d1_l{2}; case 'n' d = obj.d1_r{2}; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = d; @@ -368,8 +365,6 @@ 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 @@ -381,6 +376,7 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary case 'w' @@ -391,34 +387,32 @@ 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 % Returns the boundary sign. The right boundary is considered the positive boundary % boundary -- string function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary case {'e','n'} s = 1; case {'w','s'} s = -1; - otherwise - error('No such boundary: boundary = %s',boundary); end end % Returns borrowing constant gamma*h % boundary -- string function gamm = getBoundaryBorrowing(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + 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
--- a/+scheme/Heat2dVariable.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Heat2dVariable.m Wed Apr 10 11:00:27 2019 -0700 @@ -212,6 +212,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) if ~iscell(op) op = {op}; @@ -229,8 +230,6 @@ e = obj.e_l{2}; case 'n' e = obj.e_r{2}; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; @@ -244,8 +243,6 @@ 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 @@ -257,6 +254,7 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary case 'w' @@ -267,34 +265,32 @@ H_b = obj.H_boundary{2}; case '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) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary case {'e','n'} s = 1; case {'w','s'} s = -1; - otherwise - error('No such boundary: boundary = %s',boundary); end end % Returns borrowing constant gamma*h % boundary -- string function gamm = getBoundaryBorrowing(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + 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
--- a/+scheme/Hypsyst2d.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Hypsyst2d.m Wed Apr 10 11:00:27 2019 -0700 @@ -298,6 +298,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) if ~iscell(op) op = {op}; @@ -315,8 +316,6 @@ e = obj.e_s; case 'n' e = obj.e_n; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; end @@ -328,6 +327,7 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) e = obj.getBoundaryOperator('e', boundary); @@ -340,8 +340,6 @@ 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
--- a/+scheme/Hypsyst2dCurve.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Hypsyst2dCurve.m Wed Apr 10 11:00:27 2019 -0700 @@ -376,6 +376,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) if ~iscell(op) op = {op}; @@ -393,8 +394,6 @@ e = obj.e_s; case 'n' e = obj.e_n; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; end @@ -406,6 +405,7 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) e = obj.getBoundaryOperator('e', boundary); @@ -418,8 +418,6 @@ 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
--- a/+scheme/Laplace1d.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Laplace1d.m Wed Apr 10 11:00:27 2019 -0700 @@ -56,12 +56,13 @@ default_arg('type','neumann'); default_arg('data',0); - [e, d] = obj.getBoundaryOperator({'e', 'd'}, boundary); + e = obj.getBoundaryOperator('e', boundary); + d = obj.getBoundaryOperator('d', boundary); s = obj.getBoundarySign(boundary); switch type % Dirichlet boundary condition - case {'D','dirichlet'} + case {'D','d','dirichlet'} tuning = 1.1; tau1 = -tuning/obj.gamm; tau2 = 1; @@ -72,7 +73,7 @@ penalty = obj.a*obj.Hi*tau; % Neumann boundary condition - case {'N','neumann'} + case {'N','n','neumann'} tau = -e; closure = obj.a*obj.Hi*tau*d'; @@ -87,10 +88,12 @@ 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); + e_u = obj.getBoundaryOperator('e', boundary); + d_u = obj.getBoundaryOperator('d', boundary); s_u = obj.getBoundarySign(boundary); - [e_v, d_v] = neighbour_scheme.getBoundaryOperator({'e', 'd'}, neighbour_boundary); + e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); + d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary); s_v = neighbour_scheme.getBoundarySign(neighbour_boundary); a_u = obj.a; @@ -101,7 +104,7 @@ tuning = 1.1; - tau1 = -(a_u/gamm_u + a_v/gamm_v) * tuning; + tau1 = -1/4*(a_u/gamm_u + a_v/gamm_v) * tuning; tau2 = 1/2*a_u; sig1 = -1/2; sig2 = 0; @@ -114,51 +117,36 @@ end % Returns the boundary operator op for the boundary specified by the string boundary. - % op -- string or a cell array of strings + % op -- string % boundary -- string - function varargout = getBoundaryOperator(obj, op, boundary) + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e', 'd'}) + assertIsMember(boundary, {'l', 'r'}) - if ~iscell(op) - op = {op}; - end + o = obj.([op, '_', boundary]); + 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; + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + % Note: for 1d diffOps, the boundary quadrature is the scalar 1. + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) - 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 + H_b = 1; end % Returns the boundary sign. The right boundary is considered the positive boundary % boundary -- string function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + switch boundary case {'r'} s = 1; case {'l'} s = -1; - otherwise - error('No such boundary: boundary = %s',boundary); end end @@ -167,14 +155,4 @@ end end - - methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') - function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) - [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); - [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); - end - end -end \ No newline at end of file +end
--- a/+scheme/LaplaceCurvilinear.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/LaplaceCurvilinear.m Wed Apr 10 11:00:27 2019 -0700 @@ -278,12 +278,12 @@ default_arg('type','neumann'); default_arg('parameter', []); - e = obj.getBoundaryOperator('e', boundary); - d = obj.getBoundaryOperator('d', boundary); - H_b = obj.getBoundaryQuadrature(boundary); - s_b = obj.getBoundaryScaling(boundary); + e = obj.getBoundaryOperator('e', boundary); + d = obj.getBoundaryOperator('d', boundary); + H_b = obj.getBoundaryQuadrature(boundary); + s_b = obj.getBoundaryScaling(boundary); [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary); - m = obj.getBoundaryNumber(boundary); + m = obj.getBoundaryNumber(boundary); K = obj.K; J = obj.J; @@ -358,10 +358,10 @@ v = neighbour_scheme; % Boundary operators, u - e_u = u.getBoundaryOperator('e', boundary); - d_u = u.getBoundaryOperator('d', boundary); - gamm_u = u.getBoundaryBorrowing(boundary); - s_b_u = u.getBoundaryScaling(boundary); + e_u = u.getBoundaryOperator('e', boundary); + d_u = u.getBoundaryOperator('d', boundary); + gamm_u = u.getBoundaryBorrowing(boundary); + s_b_u = u.getBoundaryScaling(boundary); [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary); m_u = u.getBoundaryNumber(boundary); @@ -371,10 +371,10 @@ b_b_u = e_u'*u.b*e_u; % Boundary operators, v - e_v = v.getBoundaryOperator('e', neighbour_boundary); - d_v = v.getBoundaryOperator('d', neighbour_boundary); - gamm_v = v.getBoundaryBorrowing(neighbour_boundary); - s_b_v = v.getBoundaryScaling(neighbour_boundary); + e_v = v.getBoundaryOperator('e', neighbour_boundary); + d_v = v.getBoundaryOperator('d', neighbour_boundary); + gamm_v = v.getBoundaryBorrowing(neighbour_boundary); + s_b_v = v.getBoundaryScaling(neighbour_boundary); [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary); m_v = v.getBoundaryNumber(neighbour_boundary); @@ -431,16 +431,16 @@ % u denotes the solution in the own domain % v denotes the solution in the neighbour domain - e_u = obj.getBoundaryOperator('e', boundary); - d_u = obj.getBoundaryOperator('d', boundary); - H_b_u = obj.getBoundaryQuadrature(boundary); - I_u = obj.getBoundaryIndices(boundary); + e_u = obj.getBoundaryOperator('e', boundary); + d_u = obj.getBoundaryOperator('d', boundary); + H_b_u = obj.getBoundaryQuadrature(boundary); + I_u = obj.getBoundaryIndices(boundary); gamm_u = obj.getBoundaryBorrowing(boundary); - e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); - d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary); - H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary); - I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary); + e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary); + d_v = neighbour_scheme.getBoundaryOperator('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); @@ -530,6 +530,8 @@ % Returns the indices of the boundary points in the grid matrix % boundary -- string function I = getBoundaryIndices(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); switch boundary case 'w' @@ -540,14 +542,14 @@ I = ind(:,1)'; case 'n' I = ind(:,end)'; - otherwise - error('No such boundary: boundary = %s',boundary); end end % Returns borrowing constant gamma % boundary -- string function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary case {'w','e'} theta_H = obj.theta_H_u; @@ -557,8 +559,6 @@ theta_H = obj.theta_H_v; theta_M = obj.theta_M_v; theta_R = obj.theta_R_v; - otherwise - error('No such boundary: boundary = %s',boundary); end end
--- a/+scheme/Scheme.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Scheme.m Wed Apr 10 11:00:27 2019 -0700 @@ -31,20 +31,10 @@ % depending on the particular scheme implementation [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) - % TODO: op = getBoundaryOperator()?? - % makes sense to have it available through a method instead of random properties + op = getBoundaryOperator(obj, opName, boundary) + H_b= getBoundaryQuadrature(obj, boundary) % Returns the number of degrees of freedom. N = size(obj) end - - methods(Static) - % Calculates the matrcis need for the inteface coupling between - % boundary bound_u of scheme schm_u and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') - function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) - [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); - [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); - end - end end
--- a/+scheme/Schrodinger.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Schrodinger.m Wed Apr 10 11:00:27 2019 -0700 @@ -114,6 +114,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'l', 'r'}) if ~iscell(op) op = {op}; @@ -127,8 +128,6 @@ e = obj.e_l; case 'r' e = obj.e_r; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; @@ -138,24 +137,33 @@ 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 square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + % Note: for 1d diffOps, the boundary quadrature is the scalar 1. + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + + H_b = 1; + end + % Returns the boundary sign. The right boundary is considered the positive boundary % boundary -- string function s = getBoundarySign(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + switch boundary case {'r'} s = 1; case {'l'} s = -1; - otherwise - error('No such boundary: boundary = %s',boundary); end end @@ -164,14 +172,4 @@ end end - - methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') - function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) - [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); - [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); - end - end -end \ No newline at end of file +end
--- a/+scheme/Schrodinger2d.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Schrodinger2d.m Wed Apr 10 11:00:27 2019 -0700 @@ -290,6 +290,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) if ~iscell(op) op = {op}; @@ -307,8 +308,6 @@ e = obj.e_s; case 'n' e = obj.e_n; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; @@ -322,8 +321,6 @@ d = obj.d_s; case 'n' d = obj.d_n; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = d; end @@ -335,6 +332,7 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary case 'w' @@ -345,21 +343,19 @@ H_b = obj.H_boundary{2}; case '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) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary case {'e','n'} s = 1; case {'w','s'} s = -1; - otherwise - error('No such boundary: boundary = %s',boundary); end end
--- a/+scheme/TODO.txt Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,1 +0,0 @@ -% TODO: Rename package and abstract class to diffOp
--- a/+scheme/Utux.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Utux.m Wed Apr 10 11:00:27 2019 -0700 @@ -73,28 +73,24 @@ end % Returns the boundary operator op for the boundary specified by the string boundary. - % op -- string or a cell array of strings + % op -- string % boundary -- string - function varargout = getBoundaryOperator(obj, op, boundary) + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e'}) + assertIsMember(boundary, {'l', 'r'}) - if ~iscell(op) - op = {op}; - end + o = obj.([op, '_', boundary]); + 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 + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary points + % + % boundary -- string + % Note: for 1d diffOps, the boundary quadrature is the scalar 1. + function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'l', 'r'}) + + H_b = 1; end function N = size(obj) @@ -102,14 +98,4 @@ end end - - methods(Static) - % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') - function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) - [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); - [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); - end - end -end \ No newline at end of file +end
--- a/+scheme/Utux2d.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Utux2d.m Wed Apr 10 11:00:27 2019 -0700 @@ -12,6 +12,7 @@ H % Discrete norm H_x, H_y % Norms in the x and y directions Hi, Hx, Hy, Hxi, Hyi % Kroneckered norms + H_w, H_e, H_s, H_n % Boundary quadratures % Derivatives Dx, Dy @@ -59,6 +60,10 @@ Hxi = ops_x.HI; Hyi = ops_y.HI; + obj.H_w = Hy; + obj.H_e = Hy; + obj.H_s = Hx; + obj.H_n = Hx; obj.H_x = Hx; obj.H_y = Hy; obj.H = kron(Hx,Hy); @@ -272,34 +277,14 @@ end - % Returns the boundary operator op for the boundary specified by the string boundary. - % op -- string or a cell array of strings + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string % boundary -- string - function varargout = getBoundaryOperator(obj, op, boundary) - - if ~iscell(op) - op = {op}; - end + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(op, {'e'}) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) - 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 - + o = obj.([op, '_', boundary]); end % Returns square boundary quadrature matrix, of dimension @@ -307,19 +292,9 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) - 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 + H_b = obj.(['H_', boundary]); end function N = size(obj) @@ -327,14 +302,4 @@ end end - - methods(Static) - % Calculates the matrices needed for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_coupling(A,'r',B,'l') - function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) - [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); - [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); - end - end -end \ No newline at end of file +end
--- a/+scheme/Wave2d.m Fri Mar 29 14:50:50 2019 -0700 +++ b/+scheme/Wave2d.m Wed Apr 10 11:00:27 2019 -0700 @@ -200,6 +200,7 @@ % op -- string or a cell array of strings % boundary -- string function varargout = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) if ~iscell(op) op = {op}; @@ -217,8 +218,6 @@ e = obj.e_s; case 'n' e = obj.e_n; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = e; @@ -232,8 +231,6 @@ d = obj.d1_s; case 'n' d = obj.d1_n; - otherwise - error('No such boundary: boundary = %s',boundary); end varargout{i} = d; end @@ -246,6 +243,7 @@ % % boundary -- string function H_b = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) switch boundary case 'w' @@ -256,39 +254,39 @@ H_b = obj.H_x; case '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) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + 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) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + 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) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + switch boundary case 'w' Hinv = obj.Hix; @@ -298,8 +296,6 @@ Hinv = obj.Hiy; case 'n' Hinv = obj.Hiy; - otherwise - error('No such boundary: boundary = %s',boundary); end end @@ -308,14 +304,4 @@ end end - - methods(Static) - % Calculates the matrcis need for the inteface coupling between boundary bound_u of scheme schm_u - % and bound_v of scheme schm_v. - % [uu, uv, vv, vu] = inteface_couplong(A,'r',B,'l') - function [uu, uv, vv, vu] = interface_coupling(schm_u,bound_u,schm_v,bound_v) - [uu,uv] = schm_u.interface(bound_u,schm_v,bound_v); - [vv,vu] = schm_v.interface(bound_v,schm_u,bound_u); - end - end -end \ No newline at end of file +end
--- a/+scheme/Wave2dCurve.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,359 +0,0 @@ -classdef Wave2dCurve < scheme.Scheme - properties - m % Number of points in each direction, possibly a vector - h % Grid spacing - - grid - - order % Order accuracy for the approximation - - D % non-stabalized scheme operator - M % Derivative norm - c - J, Ji - a11, a12, a22 - - H % Discrete norm - Hi - H_u, H_v % Norms in the x and y directions - Hu,Hv % Kroneckerd norms. 1'*Hx*v corresponds to integration in the x dir. - Hi_u, Hi_v - Hiu, Hiv - e_w, e_e, e_s, e_n - du_w, dv_w - du_e, dv_e - du_s, dv_s - du_n, dv_n - gamm_u, gamm_v - lambda - - Dx, Dy % Physical derivatives - - x_u - x_v - y_u - y_v - end - - methods - function obj = Wave2dCurve(g ,order, c, opSet) - default_arg('opSet',@sbp.D2Variable); - default_arg('c', 1); - - warning('Use LaplaceCruveilinear instead') - - assert(isa(g, 'grid.Curvilinear')) - - m = g.size(); - m_u = m(1); - m_v = m(2); - m_tot = g.N(); - - h = g.scaling(); - h_u = h(1); - h_v = h(2); - - % Operators - ops_u = opSet(m_u, {0, 1}, order); - ops_v = opSet(m_v, {0, 1}, order); - - I_u = speye(m_u); - I_v = speye(m_v); - - D1_u = ops_u.D1; - D2_u = ops_u.D2; - H_u = ops_u.H; - Hi_u = ops_u.HI; - e_l_u = ops_u.e_l; - e_r_u = ops_u.e_r; - d1_l_u = ops_u.d1_l; - d1_r_u = ops_u.d1_r; - - D1_v = ops_v.D1; - D2_v = ops_v.D2; - H_v = ops_v.H; - Hi_v = ops_v.HI; - e_l_v = ops_v.e_l; - e_r_v = ops_v.e_r; - d1_l_v = ops_v.d1_l; - d1_r_v = ops_v.d1_r; - - Du = kr(D1_u,I_v); - Dv = kr(I_u,D1_v); - - % Metric derivatives - coords = g.points(); - x = coords(:,1); - y = coords(:,2); - - x_u = Du*x; - x_v = Dv*x; - y_u = Du*y; - y_v = Dv*y; - - J = x_u.*y_v - x_v.*y_u; - a11 = 1./J .* (x_v.^2 + y_v.^2); - a12 = -1./J .* (x_u.*x_v + y_u.*y_v); - a22 = 1./J .* (x_u.^2 + y_u.^2); - lambda = 1/2 * (a11 + a22 - sqrt((a11-a22).^2 + 4*a12.^2)); - - % Assemble full operators - L_12 = spdiags(a12, 0, m_tot, m_tot); - Duv = Du*L_12*Dv; - Dvu = Dv*L_12*Du; - - Duu = sparse(m_tot); - Dvv = sparse(m_tot); - ind = grid.funcToMatrix(g, 1:m_tot); - - for i = 1:m_v - D = D2_u(a11(ind(:,i))); - p = ind(:,i); - Duu(p,p) = D; - end - - for i = 1:m_u - D = D2_v(a22(ind(i,:))); - p = ind(i,:); - Dvv(p,p) = D; - end - - obj.H = kr(H_u,H_v); - obj.Hi = kr(Hi_u,Hi_v); - obj.Hu = kr(H_u,I_v); - obj.Hv = kr(I_u,H_v); - obj.Hiu = kr(Hi_u,I_v); - obj.Hiv = kr(I_u,Hi_v); - - obj.e_w = kr(e_l_u,I_v); - obj.e_e = kr(e_r_u,I_v); - obj.e_s = kr(I_u,e_l_v); - obj.e_n = kr(I_u,e_r_v); - obj.du_w = kr(d1_l_u,I_v); - obj.dv_w = (obj.e_w'*Dv)'; - obj.du_e = kr(d1_r_u,I_v); - obj.dv_e = (obj.e_e'*Dv)'; - obj.du_s = (obj.e_s'*Du)'; - obj.dv_s = kr(I_u,d1_l_v); - obj.du_n = (obj.e_n'*Du)'; - obj.dv_n = kr(I_u,d1_r_v); - - obj.x_u = x_u; - obj.x_v = x_v; - obj.y_u = y_u; - obj.y_v = y_v; - - obj.m = m; - obj.h = [h_u h_v]; - obj.order = order; - obj.grid = g; - - obj.c = c; - obj.J = spdiags(J, 0, m_tot, m_tot); - obj.Ji = spdiags(1./J, 0, m_tot, m_tot); - obj.a11 = a11; - obj.a12 = a12; - obj.a22 = a22; - obj.D = obj.Ji*c^2*(Duu + Duv + Dvu + Dvv); - obj.lambda = lambda; - - obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; - obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; - - obj.gamm_u = h_u*ops_u.borrowing.M.d1; - obj.gamm_v = h_v*ops_v.borrowing.M.d1; - end - - - % Closure functions return the opertors applied to the own doamin to close the boundary - % Penalty functions return the opertors to force the solution. In the case of an interface it returns the operator applied to the other doamin. - % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. - % type is a string specifying the type of boundary condition if there are several. - % data is a function returning the data that should be applied at the boundary. - % neighbour_scheme is an instance of Scheme that should be interfaced to. - % neighbour_boundary is a string specifying which boundary to interface to. - function [closure, penalty] = boundary_condition(obj, boundary, type, parameter) - default_arg('type','neumann'); - default_arg('parameter', []); - - [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv , ~, ~, ~, scale_factor] = obj.get_boundary_ops(boundary); - switch type - % Dirichlet boundary condition - case {'D','d','dirichlet'} - % v denotes the solution in the neighbour domain - tuning = 1.2; - % tuning = 20.2; - [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary); - - a_n = spdiag(coeff_n); - a_t = spdiag(coeff_t); - - F = (s * a_n * d_n' + s * a_t*d_t')'; - - u = obj; - - b1 = gamm*u.lambda./u.a11.^2; - b2 = gamm*u.lambda./u.a22.^2; - - tau = -1./b1 - 1./b2; - tau = tuning * spdiag(tau); - sig1 = 1; - - penalty_parameter_1 = halfnorm_inv_n*(tau + sig1*halfnorm_inv_t*F*e'*halfnorm_t)*e; - - closure = obj.Ji*obj.c^2 * penalty_parameter_1*e'; - penalty = -obj.Ji*obj.c^2 * penalty_parameter_1; - - - % Neumann boundary condition - case {'N','n','neumann'} - c = obj.c; - - a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); - a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); - d = (a_n * d_n' + a_t*d_t')'; - - tau1 = -s; - tau2 = 0; - tau = c.^2 * obj.Ji*(tau1*e + tau2*d); - - closure = halfnorm_inv*tau*d'; - penalty = -halfnorm_inv*tau; - - % Characteristic boundary condition - case {'characteristic', 'char', 'c'} - default_arg('parameter', 1); - beta = parameter; - c = obj.c; - - a_n = spdiags(coeff_n,0,length(coeff_n),length(coeff_n)); - a_t = spdiags(coeff_t,0,length(coeff_t),length(coeff_t)); - d = s*(a_n * d_n' + a_t*d_t')'; % outward facing normal derivative - - tau = -c.^2 * 1/beta*obj.Ji*e; - - warning('is this right?! /c?') - closure{1} = halfnorm_inv*tau/c*spdiag(scale_factor)*e'; - closure{2} = halfnorm_inv*tau*beta*d'; - penalty = -halfnorm_inv*tau; - - % Unknown, boundary condition - otherwise - error('No such boundary condition: type = %s',type); - end - end - - 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 - tuning = 1.2; - % tuning = 20.2; - [e_u, d_n_u, d_t_u, coeff_n_u, coeff_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary); - [e_v, d_n_v, d_t_v, coeff_n_v, coeff_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); - - a_n_u = spdiag(coeff_n_u); - a_t_u = spdiag(coeff_t_u); - a_n_v = spdiag(coeff_n_v); - a_t_v = spdiag(coeff_t_v); - - F_u = (s_u * a_n_u * d_n_u' + s_u * a_t_u*d_t_u')'; - F_v = (s_v * a_n_v * d_n_v' + s_v * a_t_v*d_t_v')'; - - u = obj; - v = neighbour_scheme; - - b1_u = gamm_u*u.lambda(I_u)./u.a11(I_u).^2; - b2_u = gamm_u*u.lambda(I_u)./u.a22(I_u).^2; - b1_v = gamm_v*v.lambda(I_v)./v.a11(I_v).^2; - b2_v = gamm_v*v.lambda(I_v)./v.a22(I_v).^2; - - tau = -1./(4*b1_u) -1./(4*b1_v) -1./(4*b2_u) -1./(4*b2_v); - tau = tuning * spdiag(tau); - sig1 = 1/2; - sig2 = -1/2; - - penalty_parameter_1 = halfnorm_inv_u_n*(e_u*tau + sig1*halfnorm_inv_u_t*F_u*e_u'*halfnorm_u_t*e_u); - penalty_parameter_2 = halfnorm_inv_u_n * sig2 * e_u; - - - closure = obj.Ji*obj.c^2 * ( penalty_parameter_1*e_u' + penalty_parameter_2*F_u'); - penalty = obj.Ji*obj.c^2 * (-penalty_parameter_1*e_v' + penalty_parameter_2*F_v'); - end - - % Ruturns the boundary ops and sign for the boundary specified by the string boundary. - % The right boundary is considered the positive boundary - % - % I -- the indecies of the boundary points in the grid matrix - function [e, d_n, d_t, coeff_n, coeff_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary) - - % gridMatrix = zeros(obj.m(2),obj.m(1)); - % gridMatrix(:) = 1:numel(gridMatrix); - - ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m)); - - switch boundary - case 'w' - e = obj.e_w; - d_n = obj.du_w; - d_t = obj.dv_w; - s = -1; - - I = ind(1,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); - case 'e' - e = obj.e_e; - d_n = obj.du_e; - d_t = obj.dv_e; - s = 1; - - I = ind(end,:); - coeff_n = obj.a11(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2); - case 's' - e = obj.e_s; - d_n = obj.dv_s; - d_t = obj.du_s; - s = -1; - - I = ind(:,1)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); - case 'n' - e = obj.e_n; - d_n = obj.dv_n; - d_t = obj.du_n; - s = 1; - - I = ind(:,end)'; - coeff_n = obj.a22(I); - coeff_t = obj.a12(I); - scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2); - otherwise - error('No such boundary: boundary = %s',boundary); - end - - switch boundary - case {'w','e'} - halfnorm_inv_n = obj.Hiu; - halfnorm_inv_t = obj.Hiv; - halfnorm_t = obj.Hv; - gamm = obj.gamm_u; - case {'s','n'} - halfnorm_inv_n = obj.Hiv; - halfnorm_inv_t = obj.Hiu; - halfnorm_t = obj.Hu; - gamm = obj.gamm_v; - end - end - - function N = size(obj) - N = prod(obj.m); - end - - - end -end \ No newline at end of file
--- a/+scheme/error1d.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,4 +0,0 @@ -function e = error1d(discr, v1, v2) - h = discr.h; - e = sqrt(h*sum((v1-v2).^2)); -end \ No newline at end of file
--- a/+scheme/error2d.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,5 +0,0 @@ -function e = error2d(discr, v1, v2) - % If v1 and v2 are more complex types, something like grid functions... Then we may use .getVectorFrom here! - h = discr.h; - e = sqrt(h.^2*sum((v1-v2).^2)); -end \ No newline at end of file
--- a/+scheme/errorMax.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -function e = errorMax(~, v1, v2) - e = max(abs(v1-v2)); -end
--- a/+scheme/errorRelative.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,3 +0,0 @@ -function e = errorRelative(~,v1,v2) - e = sqrt(sum((v1-v2).^2)/sum(v2.^2)); -end \ No newline at end of file
--- a/+scheme/errorSbp.m Fri Mar 29 14:50:50 2019 -0700 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,6 +0,0 @@ -function e = errorSbp(discr, v1, v2) - % If v1 and v2 are more complex types, something like grid functions... Then we may use .getVectorFrom here! - H = discr.H; - err = v2 - v1; - e = sqrt(err'*H*err); -end \ No newline at end of file