Mercurial > repos > public > sbplib
changeset 890:c70131daaa6e feature/d1_staggered
Merge with feature/poroelastic.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 21 Nov 2018 18:29:29 -0800 |
| parents | 18e10217dca9 (current diff) 7d4f57725192 (diff) |
| children | 386ef449df51 |
| files | +multiblock/multiblockgrid.m +multiblock/stitchSchemes.m +scheme/bcSetup.m +time/ExplicitRungeKuttaDiscreteData.m diracDiscr.m diracDiscrTest.m spdiagVariable.m spdiagsVariablePeriodic.m |
| diffstat | 43 files changed, 2043 insertions(+), 483 deletions(-) [+] |
line wrap: on
line diff
--- a/+anim/setup_time_quantity_plot.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+anim/setup_time_quantity_plot.m Wed Nov 21 18:29:29 2018 -0800 @@ -16,7 +16,7 @@ if ishandle(axis_handle) % t = [t t_now]; for j = 1:length(yfun) - addpoints(plot_handles(j),t_now,yfun{j}(varargin{:})); + addpoints(plot_handles(j),t_now,full(yfun{j}(varargin{:}))); end [t,~] = getpoints(plot_handles(1));
--- a/+grid/Cartesian.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+grid/Cartesian.m Wed Nov 21 18:29:29 2018 -0800 @@ -28,7 +28,11 @@ end obj.h = []; - obj.lim = []; + + obj.lim = cell(1,obj.d); + for i = 1:obj.d + obj.lim{i} = {obj.x{i}(1), obj.x{i}(end)}; + end end % n returns the number of points in the grid function o = N(obj)
--- a/+grid/evalOn.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+grid/evalOn.m Wed Nov 21 18:29:29 2018 -0800 @@ -13,18 +13,18 @@ return end % func should now be a function_handle - assert(g.D == nargin(func),'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.') + assert(g.D == nargin(func) || nargin(func) < 0,'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.') x = num2cell(g.points(),1); - k = numberOfComponents(func); + k = numberOfComponents(func, g.D); gf = func(x{:}); gf = reorderComponents(gf, k); end % Find the number of vector components of func -function k = numberOfComponents(func) - x0 = num2cell(ones(1,nargin(func))); +function k = numberOfComponents(func, dim) + x0 = num2cell(ones(1, dim)); f0 = func(x0{:}); assert(size(f0,2) == 1, 'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector'); k = length(f0);
--- a/+multiblock/DefCurvilinear.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+multiblock/DefCurvilinear.m Wed Nov 21 18:29:29 2018 -0800 @@ -48,13 +48,14 @@ g = multiblock.Grid(grids, obj.connections, obj.boundaryGroups); end - function show(obj, label, gridLines, varargin) + function h = show(obj, label, gridLines, varargin) default_arg('label', 'name') default_arg('gridLines', false); + h = []; if isempty('label') && ~gridLines for i = 1:obj.nBlocks - obj.blockMaps{i}.show(2,2); + h = [h, obj.blockMaps{i}.show(2,2)]; end axis equal return @@ -63,7 +64,7 @@ if gridLines ms = obj.getGridSizes(varargin{:}); for i = 1:obj.nBlocks - obj.blockMaps{i}.show(ms{i}(1),ms{i}(2)); + h = [h, obj.blockMaps{i}.show(ms{i}(1),ms{i}(2))]; end end @@ -76,7 +77,7 @@ for i = 1:obj.nBlocks labels{i} = num2str(i); end - otherwise + case 'none' axis equal return end
--- a/+multiblock/DiffOp.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+multiblock/DiffOp.m Wed Nov 21 18:29:29 2018 -0800 @@ -53,7 +53,11 @@ % Build the differentiation matrix - obj.blockmatrixDiv = {g.Ns, g.Ns}; + Ns = zeros(nBlocks,1); + for i = 1:nBlocks + Ns(i) = length(obj.diffOps{i}.D); + end + obj.blockmatrixDiv = {Ns, Ns}; D = blockmatrix.zero(obj.blockmatrixDiv); for i = 1:nBlocks D{i,i} = obj.diffOps{i}.D; @@ -117,7 +121,7 @@ function ops = splitOp(obj, op) % Splits a matrix operator into a cell-matrix of matrix operators for - % each g. + % each grid. ops = sparse2cell(op, obj.NNN); end @@ -144,6 +148,49 @@ end end + % Get a boundary operator specified by opName for the given boundary/BoundaryGroup + function op = getBoundaryOperatorWrapper(obj, opName, boundary) + switch class(boundary) + case 'cell' + blockId = boundary{1}; + localOp = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{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(size(obj.D,1),0); + for i = 1:length(boundary) + op = [op, obj.getBoundaryOperatorWrapper(opName, boundary{i})]; + end + otherwise + error('Unknown boundary indentifier') + end + 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); + + return + case 'multiblock.BoundaryGroup' + N = length(boundary); + H_bm = cell(N,N); + for i = 1:N + H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i}); + end + op = blockmatrix.toMatrix(H_bm); + otherwise + error('Unknown boundary indentifier') + end + end + % Creates the closure and penalty matrix for a given boundary condition, % boundary -- the name of the boundary on the form {id,name} where % id is the number of a block and name is the name of a
--- a/+multiblock/Grid.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+multiblock/Grid.m Wed Nov 21 18:29:29 2018 -0800 @@ -132,6 +132,27 @@ end + % Pads a grid function that lives on a subgrid with + % zeros and gives it the size that mathces obj. + function gf = expandFunc(obj, gfSub, subGridId) + nComponents = length(gfSub)/obj.grids{subGridId}.N(); + nBlocks = numel(obj.grids); + + % Create sparse block matrix + gfs = cell(nBlocks,1); + for i = 1:nBlocks + N = obj.grids{i}.N()*nComponents; + gfs{i} = sparse(N, 1); + end + + % Insert local gf + gfs{subGridId} = gfSub; + + % Convert cell to vector + gf = blockmatrix.toMatrix(gfs); + + end + % Find all non interface boundaries of all blocks. % Return their grid.boundaryIdentifiers in a cell array. function bs = getBoundaryNames(obj)
--- a/+multiblock/Laplace.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+multiblock/Laplace.m Wed Nov 21 18:29:29 2018 -0800 @@ -21,7 +21,7 @@ obj.D = obj.mbDiffOp.D; obj.J = obj.jacobian(); - obj.H = obj.mbDiffOp.H * obj.jacobian(); + obj.H = obj.mbDiffOp.H; end function s = size(obj) @@ -42,6 +42,10 @@ op = getBoundaryOperator(obj.mbDiffOp, opName, boundary); end + function op = getBoundaryQuadrature(obj, boundary) + op = getBoundaryQuadrature(obj.mbDiffOp, boundary); + end + function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition [closure, penalty] = boundary_condition(obj.mbDiffOp, boundary, type); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/evalOn.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,27 @@ +% Evaluate different function handle for each block in a multiblock.Grid +% Function handles may optionaly take a time argument +% f -- cell array of function handles +% f{i} = f_i(t,x,y,...) +% t -- optional time point. If not specified, it is assumed that the functions take only spatial arguments. +function gf = evalOn(g, f, t) + assertType(g, 'multiblock.Grid'); + assertType(f, 'cell'); + + default_arg('t', []); + + grids = g.grids; + nBlocks = length(grids); + gf = cell(nBlocks, 1); + + if isempty(t) + for i = 1:nBlocks + gf{i} = grid.evalOn(grids{i}, f{i}); + end + else + for i = 1:nBlocks + gf{i} = grid.evalOn(grids{i}, @(varargin)f{i}(t,varargin{:})); + end + end + + gf = blockmatrix.toMatrix(gf); +end
--- a/+multiblock/multiblockgrid.m Sun Nov 04 12:36:30 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,70 +0,0 @@ -% Creates a multi block square grid with defined boundary conditions. -% x,y defines the grid lines. Rember to think of the indexing as a matrix. Order matters! -% bc is a struct defining the boundary conditions on each side of the block. -% bc.w = {'dn',[function or value]} -function [block,conn,bound,ms] = multiblockgrid(x,y,mx,my,bc) - n = length(y)-1; % number of blocks in the y direction. - m = length(x)-1; % number of blocks in the x direction. - N = n*m; % number of blocks - - if ~issorted(x) - error('The elements of x seem to be in the wrong order'); - end - if ~issorted(flip(y)) - error('The elements of y seem to be in the wrong order'); - end - % y = sort(y,'descend'); - - % Dimensions of blocks and number of points - block = cell(n,m); - for i = 1:n - for j = 1:m - block{i,j} = { - {x(j),x(j+1)}, {y(i+1),y(i)}; - }; - - ms{i,j} = [mx(i),my(j)]; - end - end - - % Interface couplings - conn = cell(N,N); - for i = 1:n - for j = 1:m - I = flat_index(n,i,j); - if i < n - J = flat_index(n,i+1,j); - conn{I,J} = {'s','n'}; - end - - if j < m - J = flat_index(n,i,j+1); - conn{I,J} = {'e','w'}; - end - end - end - - - % Boundary conditions - bound = cell(n,m); - for i = 1:n - if isfield(bc,'w') - bound{i,1}.w = bc.w; - end - - if isfield(bc,'e') - bound{i,n}.e = bc.e; - end - end - - for j = 1:m - if isfield(bc,'n') - bound{1,j}.n = bc.n; - end - - if isfield(bc,'s') - bound{m,j}.s = bc.s; - end - end -end -
--- a/+multiblock/stitchSchemes.m Sun Nov 04 12:36:30 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,92 +0,0 @@ -% Stitch schemes together given connection matrix and BC vector. -% schmHand - function_handle to a Scheme constructor -% order - order of accuracy -% schmParam - cell array of extra parameters sent to each Scheme stored as cell arrays -% blocks - block definitions, On whatever form the scheme expects. -% ms - grid points in each direction for each block. Ex {[10,10], [10, 20]} -% conn - connection matrix -% bound - boundary condition vector, array of structs with fields w,e,s,n -% each field with a parameter array that is sent to schm.boundary_condition -% -% Output parameters are cell arrays and cell matrices. -% -% Ex: [schms, D, H] = stitchSchemes(schmHand, order, schmParam, blocks, ms, conn, bound) -function [schms, D, H] = stitchSchemes(schmHand, order, schmParam, grids, conn, bound) - default_arg('schmParam',[]); - - n_blocks = numel(grids); - - % Creating Schemes - for i = 1:n_blocks - if isempty(schmParam); - schms{i} = schmHand(grids{i},order,[]); - elseif ~iscell(schmParam) - param = schmParam(i); - schms{i} = schmHand(grids{i},order,param); - else - param = schmParam{i}; - if iscell(param) - schms{i} = schmHand(grids{i},order,param{:}); - else - schms{i} = schmHand(grids{i},order,param); - end - end - - % class(schmParam) - % class(ms) - % class(blocks) - % class(schmParam{i}) - % class(ms) - - - end - - - % Total norm - H = cell(n_blocks,n_blocks); - for i = 1:n_blocks - H{i,i} = schms{i}.H; - end - - %% Total system matrix - - % Differentiation terms - D = cell(n_blocks,n_blocks); - for i = 1:n_blocks - D{i,i} = schms{i}.D; - end - - % Boundary penalty terms - for i = 1:n_blocks - if ~isstruct(bound{i}) - continue - end - - fn = fieldnames(bound{i}); - for j = 1:length(fn); - bc = bound{i}.(fn{j}); - if isempty(bc) - continue - end - - [closure, ~] = schms{i}.boundary_condition(fn{j},bc{:}); - D{i,i} = D{i,i}+closure; - end - end - - % Interface penalty terms - for i = 1:n_blocks - for j = 1:n_blocks - intf = conn{i,j}; - if isempty(intf) - continue - end - - [uu,uv,vv,vu] = schms{i}.interface_coupling(schms{i},intf{1},schms{j},intf{2}); - D{i,i} = D{i,i} + uu; - D{i,j} = uv; - D{j,j} = D{j,j} + vv; - D{j,i} = vu; - end - end -end \ No newline at end of file
--- a/+parametrization/Curve.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+parametrization/Curve.m Wed Nov 21 18:29:29 2018 -0800 @@ -181,8 +181,8 @@ end function D = mirror(C, a, b) - assert_size(a,[2,1]); - assert_size(b,[2,1]); + assertSize(a,[2,1]); + assertSize(b,[2,1]); g = C.g; gp = C.gp; @@ -219,8 +219,8 @@ end function D = rotate(C,a,rad) - assert_size(a, [2,1]); - assert_size(rad, [1,1]); + assertSize(a, [2,1]); + assertSize(rad, [1,1]); g = C.g; gp = C.gp;
--- a/+parametrization/Ti.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+parametrization/Ti.m Wed Nov 21 18:29:29 2018 -0800 @@ -129,13 +129,13 @@ S = obj.S; if(nu>2 || nv>2) - h_grid = obj.plot(nu,nv); - set(h_grid,'Color',[0 0.4470 0.7410]); + h.grid = obj.plot(nu,nv); + set(h.grid,'Color',[0 0.4470 0.7410]); end - h_bord = obj.plot(2,2); - set(h_bord,'Color',[0.8500 0.3250 0.0980]); - set(h_bord,'LineWidth',2); + h.border = obj.plot(2,2); + set(h.border,'Color',[0.8500 0.3250 0.0980]); + set(h.border,'LineWidth',2); end
--- a/+sbp/+implementations/d2_variable_2.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/+implementations/d2_variable_2.m Wed Nov 21 18:29:29 2018 -0800 @@ -27,7 +27,7 @@ diags = -1:1; stencil = [-1/2 0 1/2]; D1 = stripeMatrix(stencil, diags, m); - + D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; D1(m,m-1)=-1;D1(m,m)=1; D1=D1/h; @@ -40,7 +40,7 @@ scheme_radius = (scheme_width-1)/2; r = (1+scheme_radius):(m-scheme_radius); - function D2 = D2_fun(c) + function [D2, B] = D2_fun(c) Mm1 = -c(r-1)/2 - c(r)/2; M0 = c(r-1)/2 + c(r) + c(r+1)/2; @@ -54,6 +54,8 @@ M=M/h; D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); + B = HI*M; end D2 = @D2_fun; + end \ No newline at end of file
--- a/+sbp/+implementations/d2_variable_4.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/+implementations/d2_variable_4.m Wed Nov 21 18:29:29 2018 -0800 @@ -49,7 +49,7 @@ N = m; - function D2 = D2_fun(c) + function [D2, B] = D2_fun(c) M = 78+(N-12)*5; %h = 1/(N-1); @@ -131,6 +131,8 @@ cols(40+(i-7)*5:44+(i-7)*5) = [i-2;i-1;i;i+1;i+2]; end D2 = sparse(rows,cols,D2); + + B = HI*( c(end)*e_r*d1_r' - c(1)*e_l*d1_l') - D2; end D2 = @D2_fun; end \ No newline at end of file
--- a/+sbp/+implementations/d2_variable_periodic_2.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/+implementations/d2_variable_periodic_2.m Wed Nov 21 18:29:29 2018 -0800 @@ -27,7 +27,7 @@ scheme_width = 3; scheme_radius = (scheme_width-1)/2; - + r = 1:m; offset = scheme_width; r = r + offset; @@ -41,7 +41,7 @@ vals = [Mm1,M0,Mp1]; diags = -scheme_radius : scheme_radius; - M = spdiagsVariablePeriodic(vals,diags); + M = spdiagsPeriodic(vals,diags); M=M/h; D2=HI*(-M );
--- a/+sbp/+implementations/d2_variable_periodic_4.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/+implementations/d2_variable_periodic_4.m Wed Nov 21 18:29:29 2018 -0800 @@ -30,7 +30,7 @@ scheme_width = 5; scheme_radius = (scheme_width-1)/2; - + r = 1:m; offset = scheme_width; r = r + offset; @@ -47,7 +47,7 @@ vals = -[Mm2,Mm1,M0,Mp1,Mp2]; diags = -scheme_radius : scheme_radius; - M = spdiagsVariablePeriodic(vals,diags); + M = spdiagsPeriodic(vals,diags); M=M/h; D2=HI*(-M );
--- a/+sbp/+implementations/d2_variable_periodic_6.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/+implementations/d2_variable_periodic_6.m Wed Nov 21 18:29:29 2018 -0800 @@ -47,12 +47,12 @@ vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3]; diags = -scheme_radius : scheme_radius; - M = spdiagsVariablePeriodic(vals,diags); + M = spdiagsPeriodic(vals,diags); M=M/h; D2=HI*(-M ); end D2 = @D2_fun; - + end
--- a/+sbp/+implementations/d4_variable_6.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/+implementations/d4_variable_6.m Wed Nov 21 18:29:29 2018 -0800 @@ -85,7 +85,7 @@ scheme_radius = (scheme_width-1)/2; r = (1+scheme_radius):(m-scheme_radius); - function D2 = D2_fun(c) + function [D2, B] = D2_fun(c) Mm3 = c(r-2)/0.40e2 + c(r-1)/0.40e2 - 0.11e2/0.360e3 * c(r-3) - 0.11e2/0.360e3 * c(r); Mm2 = c(r-3)/0.20e2 - 0.3e1/0.10e2 * c(r-1) + c(r+1)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r-2); @@ -128,6 +128,7 @@ M=M/h; D2 = HI*(-M - c(1)*e_l*d1_l' + c(m)*e_r*d1_r'); + B = HI*M; end D2 = @D2_fun;
--- a/+sbp/InterpAWW.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+sbp/InterpAWW.m Wed Nov 21 18:29:29 2018 -0800 @@ -23,6 +23,7 @@ % accOp : String, 'C2F' or 'F2C'. Specifies which of the operators % should have higher accuracy. function obj = InterpAWW(m_C,m_F,order_C,order_F,accOp) + assertIsMember(accOp, {'C2F','F2C'}); ratio = (m_F-1)/(m_C-1); h_C = 1;
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Elastic2dCurvilinear.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,621 @@ +classdef Elastic2dCurvilinear < scheme.Scheme + +% 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 +% +% 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 +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrices for varible coefficients + LAMBDA % Variable coefficient, related to dilation + MU % Shear modulus, variable coefficient + RHO, RHOi % Density, variable + + % Metric coefficients + b % Cell matrix of size dim x dim + J, Ji + beta % Cell array of scale factors + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_lambda + D2_mu + + % Traction operators used for BC + 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 + 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 + RHOi_kron + Ji_kron, J_kron + Hi_kron, H_kron + end + + methods + + function obj = Elastic2dCurvilinear(g ,order, lambda_fun, mu_fun, rho_fun, 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); + dim = 2; + + lambda = grid.evalOn(g, lambda_fun); + mu = grid.evalOn(g, mu_fun); + rho = grid.evalOn(g, rho_fun); + m = g.size(); + obj.m = m; + m_tot = g.N(); + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), {0, 1}, order); + end + + % 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; + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + %====== Assemble full operators ======== + + % Variable coefficients + LAMBDA = spdiag(lambda); + obj.LAMBDA = LAMBDA; + MU = spdiag(mu); + obj.MU = MU; + RHO = spdiag(rho); + obj.RHO = RHO; + obj.RHOi = inv(RHO); + + % Allocate + obj.D1 = cell(dim,1); + obj.D2_lambda = cell(dim,dim,dim); + obj.D2_mu = cell(dim,dim,dim); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % -- Metric coefficients ---- + coords = g.points(); + x = coords(:,1); + y = coords(:,2); + + % Use non-periodic difference operators for metric even if opSet is periodic. + xmax = max(ops{1}.x); + ymax = max(ops{2}.x); + 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); + + x_xi = D1Metric{1}*x; + x_eta = D1Metric{2}*x; + y_xi = D1Metric{1}*y; + y_eta = D1Metric{2}*y; + + J = x_xi.*y_eta - x_eta.*y_xi; + + b = cell(dim,dim); + b{1,1} = y_eta./J; + b{1,2} = -x_eta./J; + b{2,1} = -y_xi./J; + b{2,2} = x_xi./J; + + % Scale factors for boundary integrals + beta = cell(dim,1); + beta{1} = sqrt(x_eta.^2 + y_eta.^2); + beta{2} = sqrt(x_xi.^2 + y_xi.^2); + + J = spdiag(J); + Ji = inv(J); + for i = 1:dim + beta{i} = spdiag(beta{i}); + for j = 1:dim + b{i,j} = spdiag(b{i,j}); + end + end + obj.J = J; + obj.Ji = Ji; + obj.b = b; + obj.beta = beta; + %---------------------------- + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 + for i = 1:dim + for j = 1:dim + for k = 1:dim + obj.D2_lambda{i,j,k} = sparse(m_tot); + obj.D2_mu{i,j,k} = sparse(m_tot); + end + end + end + ind = grid.funcToMatrix(g, 1:m_tot); + + % x-dir + for i = 1:dim + for j = 1:dim + for k = 1 + + coeff_lambda = J*b{i,k}*b{j,k}*lambda; + coeff_mu = J*b{j,k}*b{i,k}*mu; + + for col = 1:m(2) + D_lambda = D2{1}(coeff_lambda(ind(:,col))); + D_mu = D2{1}(coeff_mu(ind(:,col))); + + p = ind(:,col); + obj.D2_lambda{i,j,k}(p,p) = D_lambda; + obj.D2_mu{i,j,k}(p,p) = D_mu; + end + + end + end + end + + % y-dir + for i = 1:dim + for j = 1:dim + for k = 2 + + coeff_lambda = J*b{i,k}*b{j,k}*lambda; + coeff_mu = J*b{j,k}*b{i,k}*mu; + + for row = 1:m(1) + D_lambda = D2{2}(coeff_lambda(ind(row,:))); + D_mu = D2{2}(coeff_mu(ind(row,:))); + + p = ind(row,:); + obj.D2_lambda{i,j,k}(p,p) = D_lambda; + obj.D2_mu{i,j,k}(p,p) = D_mu; + end + + end + end + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary_l = cell(dim,1); + obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2}; + obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1}; + obj.H_boundary_r = cell(dim,1); + obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2}; + obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1}; + + % E{i}^T picks out component i. + E = cell(dim,1); + I = speye(m_tot,m_tot); + for i = 1:dim + e = sparse(dim,1); + e(i) = 1; + E{i} = kron(I,e); + end + obj.E = E; + + % Differentiation matrix D (without SAT) + D2_lambda = obj.D2_lambda; + D2_mu = obj.D2_mu; + D1 = obj.D1; + D = sparse(dim*m_tot,dim*m_tot); + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_lambda{i,j,k}*E{j}' + ... + db(k,l)*D1{k}*J*b{i,k}*b{j,l}*LAMBDA*D1{l}*E{j}' ... + ); + + D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{i,j,k}*E{j}' + ... + db(k,l)*D1{k}*J*b{j,k}*b{i,l}*MU*D1{l}*E{j}' ... + ); + + D = D + E{i}*Ji*inv(RHO)*( d(k,l)*D2_mu{j,j,k}*E{i}' + ... + db(k,l)*D1{k}*J*b{j,k}*b{j,l}*MU*D1{l}*E{i}' ... + ); + + end + end + end + end + obj.D = D; + %=========================================% + + % Numerical traction operators for BC. + % Because d1 =/= e0^T*D1, the numerical tractions are different + % at every boundary. + T_l = cell(dim,1); + T_r = cell(dim,1); + tau_l = cell(dim,1); + tau_r = cell(dim,1); + % tau^{j}_i = sum_k T^{j}_{ik} u_k + + d1_l = obj.d1_l; + d1_r = obj.d1_r; + e_l = obj.e_l; + e_r = obj.e_r; + + % Loop over boundaries + for j = 1:dim + T_l{j} = cell(dim,dim); + T_r{j} = cell(dim,dim); + tau_l{j} = cell(dim,1); + tau_r{j} = cell(dim,1); + + % 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); + + % Loop over components that T_{ik}^{(j)} acts on + for k = 1:dim + + T_l{j}{i,k} = sparse(m_tot,m_tot); + T_r{j}{i,k} = sparse(m_tot,m_tot); + + for m = 1:dim + for l = 1:dim + 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} + ... + 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}); + end + 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}; + + 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}'; + end + + end + end + obj.T_l = T_l; + obj.T_r = T_r; + obj.tau_l = tau_l; + obj.tau_r = tau_r; + + % Kroneckered norms and coefficients + I_dim = speye(dim); + obj.RHOi_kron = kron(obj.RHOi, I_dim); + obj.Ji_kron = kron(obj.Ji, I_dim); + obj.Hi_kron = kron(obj.Hi, I_dim); + obj.H_kron = kron(obj.H, I_dim); + obj.J_kron = kron(obj.J, I_dim); + + % Misc. + obj.h = g.scaling(); + obj.order = order; + obj.grid = g; + obj.dim = dim; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators 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'. + % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition + % on the first component. + % 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, bc, tuning) + default_arg('tuning', 1.2); + + assert( iscell(bc), 'The BC type must be a 2x1 cell array' ); + comp = bc{1}; + type = bc{2}; + + % 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 = obj.E; + Hi = obj.Hi; + LAMBDA = obj.LAMBDA; + MU = obj.MU; + RHOi = obj.RHOi; + Ji = obj.Ji; + + dim = obj.dim; + m_tot = obj.grid.N(); + + % Preallocate + closure = sparse(dim*m_tot, dim*m_tot); + penalty = sparse(dim*m_tot, m_tot/obj.m(j)); + + % Loop over components that we (potentially) have different BC on + k = comp; + switch type + + % 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 ); + + % 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}' ); + penalty = penalty - E{i}*RHOi*Hi*Ji*B'*e*H_gamma; + 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} ); + penalty = penalty + E{k}*RHOi*Ji*Hi*e*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + % Operators without subscripts are from the own domain. + error('Not implemented'); + 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); + + % 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; + + 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; + %------------------------- + + % Borrowing constants + phi_u = obj.phi{j}; + h_u = obj.h(j); + h11_u = obj.H11{j}*h_u; + gamma_u = obj.gamma{j}; + + phi_v = neighbour_scheme.phi{j_v}; + h_v = neighbour_scheme.h(j_v); + h11_v = neighbour_scheme.H11{j_v}*h_v; + 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) + th1 = h11/(2*dim); + th2 = h11*phi/2; + th3 = h*gamma; + a1 = ( (th1 + th2)*th3*lambda + 4*th1*th2*mu ) / (2*th1*th2*th3); + a2 = ( 16*(th1 + th2)*lambda*mu ) / (th1*th2*th3); + alpha_ii = a1 + sqrt(a2 + a1^2); + + alpha_ij = mu*(2/h11 + 1/(phi*h11)); + 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); + sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4; + sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4; + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij); + + % 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*sigma(i,j)*H_gamma*e'*E{i}'; + penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}'; + + 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}; + + % 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 + 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) + + 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 + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + 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 + function [varargout] = get_boundary_operator(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); + end + + if ~iscell(op) + op = {op}; + end + + for i = 1:length(op) + switch op{i} + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.d1_r{j}; + end + case 'H' + 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'} + varargout{i} = obj.H_boundary_r{j}; + end + case 'T' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.T_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.T_r{j}; + end + case 'tau' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.tau_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.tau_r{j}; + end + otherwise + error(['No such operator: operator = ' op{i}]); + end + end + + end + + function N = size(obj) + N = obj.dim*prod(obj.m); + end + end +end
--- a/+scheme/Elastic2dVariable.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+scheme/Elastic2dVariable.m Wed Nov 21 18:29:29 2018 -0800 @@ -1,7 +1,7 @@ classdef Elastic2dVariable < scheme.Scheme % Discretizes the elastic wave 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 % opSet should be cell array of opSets, one per dimension. This % is useful if we have periodic BC in one direction. @@ -30,40 +30,60 @@ 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 d1_l, d1_r % Normal derivatives at the boundary E % E{i}^T picks out component i - + H_boundary % Boundary inner products % 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(); h = g.scaling(); lim = g.lim; + if isempty(lim) + x = g.x; + lim = cell(length(x),1); + for i = 1:length(x) + lim{i} = {min(x{i}), max(x{i})}; + end + end % 1D operators ops = cell(dim,1); @@ -73,10 +93,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; + 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); @@ -164,6 +183,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); @@ -194,7 +214,7 @@ end end obj.D = D; - %=========================================% + %=========================================%' % Numerical traction operators for BC. % Because d1 =/= e0^T*D1, the numerical tractions are different @@ -223,14 +243,14 @@ tau_l{j}{i} = sparse(m_tot,dim*m_tot); tau_r{j}{i} = sparse(m_tot,dim*m_tot); for k = 1:dim - T_l{j}{i,k} = ... + 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(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}'; - T_r{j}{i,k} = ... + 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(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}'; tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; @@ -256,41 +276,41 @@ obj.grid = g; obj.dim = dim; + % Used for adjoint optimization + obj.B = cell(1,dim); + for i = 1:dim + obj.B{i} = zeros(m(i),m(i),m(i)); + for k = 1:m(i) + c = sparse(m(i),1); + c(k) = 1; + [~, obj.B{i}(:,:,k)] = ops{i}.D2(c); + end + end + end % Closure functions return the operators applied to the own domain to close the boundary % Penalty functions return the operators 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 cell array of strings specifying the type of boundary condition for each component. + % bc is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition + % on the first component. % 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',{'free','free'}); - default_arg('parameter', []); + function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) + default_arg('tuning', 1.2); + + assert( iscell(bc), 'The BC type must be a 2x1 cell array' ); + comp = bc{1}; + type = bc{2}; % 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; - tau = obj.tau_r{j}; - T = obj.T_r{j}; - case -1 - e = obj.e_l; - d = obj.d1_l; - tau = obj.tau_l{j}; - T = obj.T_l{j}; - end + j = obj.get_boundary_number(boundary); + [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); E = obj.E; Hi = obj.Hi; - H_gamma = obj.H_boundary{j}; LAMBDA = obj.LAMBDA; MU = obj.MU; RHOi = obj.RHOi; @@ -298,66 +318,126 @@ dim = obj.dim; m_tot = obj.grid.N(); - RHOi_kron = obj.RHOi_kron; - Hi_kron = obj.Hi_kron; - % Preallocate closure = sparse(dim*m_tot, dim*m_tot); - penalty = cell(dim,1); - for k = 1:dim - penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j)); - end + penalty = sparse(dim*m_tot, m_tot/obj.m(j)); - % Loop over components that we (potentially) have different BC on - for k = 1:dim - switch type{k} + k = comp; + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} - % Dirichlet boundary condition - case {'D','d','dirichlet','Dirichlet'} + theta_R = obj.theta_R{j}; + theta_H = obj.theta_H{j}; + theta_M = obj.theta_M{j}; - tuning = 1.2; - 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); + 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 = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... - + d(i,j)* a_mu_i*MU ... - + db(i,j)*a_mu_ij*MU ); + 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 ); - % 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*B'*e{j}*H_gamma*(e{j}'*E{k}' ); - penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma; - end + % 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*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{j}*H_gamma* (e{j}'*tau{k} ); - penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma; + % Free boundary condition + case {'F','f','Free','free','traction','Traction','t','T'} + closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); + penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; - % Unknown boundary condition - otherwise - error('No such boundary condition: type = %s',type); - 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) % 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; - % tuning = 20.2; - error('Interface not implemented'); + + % 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); + + % 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; + + 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; + %------------------------- + + % Borrowing constants + theta_R_u = obj.theta_R{j}; + theta_H_u = obj.theta_H{j}; + theta_M_u = obj.theta_M{j}; + + theta_R_v = neighbour_scheme.theta_R{j_v}; + theta_H_v = neighbour_scheme.theta_H{j_v}; + theta_M_v = neighbour_scheme.theta_M{j_v}; + + function [alpha_ii, alpha_ij] = computeAlpha(th_R, th_H, th_M, lambda, mu) + alpha_ii = dim*lambda/(4*th_H) + lambda/(4*th_R) + mu/(2*th_M); + alpha_ij = mu/(2*th_H) + mu/(4*th_R); + end + + [alpha_ii_u, alpha_ij_u] = computeAlpha(theta_R_u, theta_H_u, theta_M_u, lambda_u, mu_u); + [alpha_ii_v, alpha_ij_v] = computeAlpha(theta_R_v, theta_H_v, theta_M_v, lambda_v, mu_v); + sigma_ii = tuning*(alpha_ii_u + alpha_ii_v); + sigma_ij = tuning*(alpha_ij_u + alpha_ij_v); + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij); + + % 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*sigma(i,j)*H_gamma*e'*E{i}'; + penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}'; + + 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}; + + % 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 end % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. @@ -380,8 +460,9 @@ 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 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) switch boundary case {'w','W','west','West', 'e', 'E', 'east', 'East'} @@ -392,29 +473,73 @@ error('No such boundary: boundary = %s',boundary); 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 - otherwise - error(['No such operator: operatr = ' op]); + if ~iscell(op) + op = {op}; + end + + for i = 1:length(op) + switch op{i} + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = 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}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.T_r{j}; + end + case 'tau' + switch boundary + case {'w','W','west','West','s','S','south','South'} + varargout{i} = obj.tau_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{i} = obj.tau_r{j}; + end + case 'alpha' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + tuning = 1.2; + LAMBDA = obj.LAMBDA; + MU = obj.MU; + + phi = obj.phi{j}; + h = obj.h(j); + h11 = obj.H11{j}*h; + gamma = obj.gamma{j}; + dim = obj.dim; + + 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,k) d(i,k)*tuning*( d(i,j)* a_lambda*LAMBDA ... + + d(i,j)* a_mu_i*MU ... + + db(i,j)*a_mu_ij*MU ); + varargout{i} = alpha; + otherwise + error(['No such operator: operator = ' op{i}]); + end end end function N = size(obj) - N = prod(obj.m); + N = obj.dim*prod(obj.m); end end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Heat2dCurvilinear.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,385 @@ +classdef Heat2dCurvilinear < scheme.Scheme + +% Discretizes the Laplacian with variable coefficent, curvilinear, +% in the Heat equation way (i.e., the discretization matrix is not necessarily +% symmetric) +% 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. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrix for variable coefficients + KAPPA % Variable coefficient + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_kappa + + H, Hi % Inner products + 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 + + % Metric coefficients + b % Cell matrix of size dim x dim + J, Ji + beta % Cell array of scale factors + + % Numerical boundary flux operators + flux_l, flux_r + + end + + methods + + function obj = Heat2dCurvilinear(g ,order, kappa_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('kappa_fun', @(x,y) 0*x+1); + dim = 2; + + kappa = grid.evalOn(g, kappa_fun); + m = g.size(); + m_tot = g.N(); + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), {0, 1}, order); + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + %====== Assemble full operators ======== + KAPPA = spdiag(kappa); + obj.KAPPA = KAPPA; + + % Allocate + obj.D1 = cell(dim,1); + obj.D2_kappa = cell(dim,1); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % -- Metric coefficients ---- + coords = g.points(); + x = coords(:,1); + y = coords(:,2); + + % Use non-periodic difference operators for metric even if opSet is periodic. + xmax = max(ops{1}.x); + ymax = max(ops{2}.x); + 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); + + x_xi = D1Metric{1}*x; + x_eta = D1Metric{2}*x; + y_xi = D1Metric{1}*y; + y_eta = D1Metric{2}*y; + + J = x_xi.*y_eta - x_eta.*y_xi; + + b = cell(dim,dim); + b{1,1} = y_eta./J; + b{1,2} = -x_eta./J; + b{2,1} = -y_xi./J; + b{2,2} = x_xi./J; + + % Scale factors for boundary integrals + beta = cell(dim,1); + beta{1} = sqrt(x_eta.^2 + y_eta.^2); + beta{2} = sqrt(x_xi.^2 + y_xi.^2); + + J = spdiag(J); + Ji = inv(J); + for i = 1:dim + beta{i} = spdiag(beta{i}); + for j = 1:dim + b{i,j} = spdiag(b{i,j}); + end + end + obj.J = J; + obj.Ji = Ji; + obj.b = b; + obj.beta = beta; + %---------------------------- + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 coefficients + kappa_coeff = cell(dim,dim); + for j = 1:dim + 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; + end + end + ind = grid.funcToMatrix(g, 1:m_tot); + + % x-dir + j = 1; + for col = 1:m(2) + D_kappa = D2{1}(kappa_coeff{j}(ind(:,col))); + + p = ind(:,col); + obj.D2_kappa{j}(p,p) = D_kappa; + end + + % y-dir + j = 2; + for row = 1:m(1) + D_kappa = D2{2}(kappa_coeff{j}(ind(row,:))); + + p = ind(row,:); + obj.D2_kappa{j}(p,p) = D_kappa; + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary_l = cell(dim,1); + obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2}; + obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1}; + obj.H_boundary_r = cell(dim,1); + obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2}; + obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1}; + + %=== Differentiation matrix D (without SAT) === + D2_kappa = obj.D2_kappa; + D1 = obj.D1; + D = sparse(m_tot,m_tot); + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + + % 2nd derivatives + for j = 1:dim + D = D + Ji*D2_kappa{j}; + end + + % Mixed terms + for i = 1:dim + for j = 1:dim + for k = 1:dim + D = D + db(i,j)*Ji*D1{j}*b{i,j}*J*KAPPA*b{i,k}*D1{k}; + end + end + end + obj.D = D; + %=========================================% + + % Normal flux operators for BC. + flux_l = cell(dim,1); + flux_r = cell(dim,1); + + d1_l = obj.d1_l; + d1_r = obj.d1_r; + e_l = obj.e_l; + e_r = obj.e_r; + + % Loop over boundaries + for j = 1:dim + flux_l{j} = sparse(m_tot,m_tot); + flux_r{j} = sparse(m_tot,m_tot); + + % Loop over dummy index + for i = 1:dim + % Loop over dummy index + for k = 1:dim + flux_l{j} = flux_l{j} ... + - beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_l{k}*d1_l{k}' + db(j,k)*D1{k} ); + + flux_r{j} = flux_r{j} ... + + beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_r{k}*d1_r{k}' + db(j,k)*D1{k} ); + end + + end + end + obj.flux_l = flux_l; + obj.flux_r = flux_r; + + % Misc. + obj.m = m; + obj.h = g.scaling(); + obj.order = order; + obj.grid = g; + obj.dim = dim; + obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1]; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators 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. + % 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, symmetric, tuning) + default_arg('type','Neumann'); + default_arg('symmetric', false); + default_arg('tuning',1.2); + + % 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{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 + + Hi = obj.Hi; + Ji = obj.Ji; + KAPPA = obj.KAPPA; + kappa_gamma = e'*KAPPA*e; + h = obj.h(j); + alpha = h*obj.alpha(j); + + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + + if ~symmetric + 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' ) ; + 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; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + 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) + + 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 + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + nj = 1; + 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) + + 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 + + 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 + otherwise + error(['No such operator: operatr = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end
--- a/+scheme/Heat2dVariable.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+scheme/Heat2dVariable.m Wed Nov 21 18:29:29 2018 -0800 @@ -28,6 +28,7 @@ H, Hi % Inner products e_l, e_r d1_l, d1_r % Normal derivatives at the boundary + alpha % Vector of borrowing constants H_boundary % Boundary inner products @@ -144,6 +145,7 @@ obj.order = order; obj.grid = g; obj.dim = dim; + obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1]; end @@ -155,9 +157,10 @@ % 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) + function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning) default_arg('type','Neumann'); - default_arg('parameter', []); + default_arg('symmetric', false); + default_arg('tuning',1.2); % j is the coordinate direction of the boundary % nj: outward unit normal component. @@ -177,18 +180,29 @@ H_gamma = obj.H_boundary{j}; KAPPA = obj.KAPPA; kappa_gamma = e{j}'*KAPPA*e{j}; + h = obj.h(j); + alpha = h*obj.alpha(j); switch type % Dirichlet boundary condition 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; + 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; + end % Free boundary condition case {'N','n','neumann','Neumann'} closure = -nj*Hi*e{j}*kappa_gamma*H_gamma*(d{j}' ); - penalty = nj*Hi*e{j}*kappa_gamma*H_gamma; + penalty = Hi*e{j}*kappa_gamma*H_gamma; + % penalty is for normal derivative and not for derivative, hence the sign. % Unknown boundary condition otherwise
--- a/+scheme/bcSetup.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+scheme/bcSetup.m Wed Nov 21 18:29:29 2018 -0800 @@ -11,49 +11,34 @@ % In the case where it only depends on time it should return the data as grid function for the boundary. % In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain. % For example in the 2D case: f(t,x,y). -function [closure, S] = bcSetup(diffOp, bc, S_sign) +function [closure, S] = bcSetup(diffOp, bcs, S_sign) default_arg('S_sign', 1); - assertType(bc, 'cell'); + assertType(bcs, 'cell'); assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1'); + verifyBcFormat(bcs, diffOp); % Setup storage arrays closure = spzeros(size(diffOp)); - gridDataPenalties = {}; - gridDataFunctions = {}; - symbolicDataPenalties = {}; - symbolicDataFunctions = {}; - symbolicDataCoords = {}; + gridData = {}; + symbolicData = {}; % Collect closures, penalties and data - for i = 1:length(bc) - assertType(bc{i}, 'struct'); - [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type); + for i = 1:length(bcs) + [localClosure, penalty] = diffOp.boundary_condition(bcs{i}.boundary, bcs{i}.type); closure = closure + localClosure; - if ~isfield(bc{i},'data') || isempty(bc{i}.data) - % Skip to next loop if there is no data + [ok, isSymbolic, data] = parseData(bcs{i}, penalty, diffOp.grid); + + if ~ok + % There was no data continue end - assertType(bc{i}.data, 'function_handle'); - % Find dimension - dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2); - - if nargin(bc{i}.data) == 1 - % Grid data - boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1]; - assert_size(bc{i}.data(0), boundarySize); % Eval for t = 0 and make sure the function returns a grid vector of the correct size. - gridDataPenalties{end+1} = penalty; - gridDataFunctions{end+1} = bc{i}.data; - elseif nargin(bc{i}.data) == 1+dim - % Symbolic data - coord = diffOp.grid.getBoundary(bc{i}.boundary); - symbolicDataPenalties{end+1} = penalty; - symbolicDataFunctions{end+1} = bc{i}.data; - symbolicDataCoords{end+1} = num2cell(coord ,1); + if isSymbolic + symbolicData{end+1} = data; else - error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bc{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i); + gridData{end+1} = data; end end @@ -61,15 +46,67 @@ O = spzeros(size(diffOp),1); function v = S_fun(t) v = O; - for i = 1:length(gridDataFunctions) - v = v + gridDataPenalties{i}*gridDataFunctions{i}(t); + for i = 1:length(gridData) + v = v + gridData{i}.penalty*gridData{i}.func(t); end - for i = 1:length(symbolicDataFunctions) - v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:}); + for i = 1:length(symbolicData) + v = v + symbolicData{i}.penalty*symbolicData{i}.func(t, symbolicData{i}.coords{:}); end v = S_sign * v; end S = @S_fun; end + +function verifyBcFormat(bcs, diffOp) + for i = 1:length(bcs) + assertType(bcs{i}, 'struct'); + assertStructFields(bcs{i}, {'type', 'boundary'}); + + if ~isfield(bcs{i}, 'data') || isempty(bcs{i}.data) + continue + end + + if ~isa(bcs{i}.data, 'function_handle') + error('bcs{%d}.data should be a function of time or a function of time and space',i); + end + + b = diffOp.grid.getBoundary(bcs{i}.boundary); + + dim = size(b,2); + + if nargin(bcs{i}.data) == 1 + % Grid data (only function of time) + assertSize(bcs{i}.data(0), 1, size(b)); + elseif nargin(bcs{i}.data) ~= 1+dim + error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bcs{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i); + end + end +end + +function [ok, isSymbolic, dataStruct] = parseData(bc, penalty, grid) + if ~isfield(bc,'data') || isempty(bc.data) + isSymbolic = []; + dataStruct = struct(); + ok = false; + return + end + ok = true; + + nArg = nargin(bc.data); + + if nArg > 1 + % Symbolic data + isSymbolic = true; + coord = grid.getBoundary(bc.boundary); + dataStruct.penalty = penalty; + dataStruct.func = bc.data; + dataStruct.coords = num2cell(coord, 1); + else + % Grid data + isSymbolic = false; + dataStruct.penalty = penalty; + dataStruct.func = bcs{i}.data; + end +end
--- a/+time/ExplicitRungeKuttaDiscreteData.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+time/ExplicitRungeKuttaDiscreteData.m Wed Nov 21 18:29:29 2018 -0800 @@ -3,7 +3,6 @@ D S % Function handle for time-dependent data data % Matrix of data vectors, one column per stage - F k t v @@ -57,7 +56,13 @@ a = obj.a; b = obj.b; c = obj.c; - s = obj.s; + s = obj.s; + end + + % Returns quadrature weights for stages in one time step + function quadWeights = getTimeStepQuadrature(obj) + [~, b] = obj.getTableau(); + quadWeights = obj.k*b; end function obj = step(obj) @@ -82,12 +87,12 @@ K(:,i) = D*U(:,i); obj.T(i) = obj.t + c(i)*dt; - % Data from continuos function and discrete time-points. + % Data from continuous function and discrete time-points. if ~isempty(S) - K(:,i) = K(:,i) + S(obj.T(i)); + K(:,i) = K(:,i) + S(obj.T(i)); end if ~isempty(data) - K(:,i) = K(:,i) + data(:,obj.n*s + i); + K(:,i) = K(:,i) + data(:,obj.n*s + i); end end @@ -102,11 +107,11 @@ methods (Static) - function k = getTimeStep(lambda) - - switch obj.order + function k = getTimeStep(lambda, order) + default_arg('order', 4); + switch order case 4 - k = rk4.get_rk4_time_step(lambda); + k = time.rk4.get_rk4_time_step(lambda); otherwise error('Time-step function not available for this order'); end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/ExplicitRungeKuttaSecondOrderDiscreteData.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,129 @@ +classdef ExplicitRungeKuttaSecondOrderDiscreteData < time.Timestepper + properties + k + t + w + m + D + E + M + C_cont % Continuous part (function handle) of forcing on first order form. + C_discr% Discrete part (matrix) of forcing on first order form. + n + order + tsImplementation % Time stepper object, RK first order form, + % which this wraps around. + end + + + methods + % Solves u_tt = Du + Eu_t + S by + % Rewriting on first order form: + % w_t = M*w + C(t) + % where + % M = [ + % 0, I; + % D, E; + % ] + % and + % C(t) = [ + % 0; + % S(t) + % ] + % D, E, should be matrices (or empty for zero) + % They can also be omitted by setting them equal to the empty matrix. + % S = S_cont + S_discr, where S_cont is a function handle + % S_discr a matrix of data vectors, one column per stage. + function obj = ExplicitRungeKuttaSecondOrderDiscreteData(D, E, S_cont, S_discr, k, t0, v0, v0t, order) + default_arg('order', 4); + default_arg('S_cont', []); + default_arg('S_discr', []); + obj.D = D; + obj.E = E; + obj.m = length(v0); + obj.n = 0; + + default_arg('D', sparse(obj.m, obj.m) ); + default_arg('E', sparse(obj.m, obj.m) ); + + obj.k = k; + obj.t = t0; + obj.w = [v0; v0t]; + + I = speye(obj.m); + O = sparse(obj.m,obj.m); + + obj.M = [ + O, I; + D, E; + ]; + + % Build C_cont + if ~isempty(S_cont) + obj.C_cont = @(t)[ + sparse(obj.m,1); + S_cont(t) + ]; + else + obj.C_cont = []; + end + + % Build C_discr + if ~isempty(S_discr) + [~, nt] = size(S_discr); + obj.C_discr = [sparse(obj.m, nt); + S_discr + ]; + else + obj.C_discr = []; + end + obj.tsImplementation = time.ExplicitRungeKuttaDiscreteData(obj.M, obj.C_cont, obj.C_discr,... + k, obj.t, obj.w, order); + end + + function [v,t,U,T,K] = getV(obj) + [w,t,U,T,K] = obj.tsImplementation.getV(); + + v = w(1:end/2); + U = U(1:end/2, :); % Stage approximations in previous time step. + K = K(1:end/2, :); % Stage rates in previous time step. + % T: Stage times in previous time step. + end + + function [vt,t,U,T,K] = getVt(obj) + [w,t,U,T,K] = obj.tsImplementation.getV(); + + vt = w(end/2+1:end); + U = U(end/2+1:end, :); % Stage approximations in previous time step. + K = K(end/2+1:end, :); % Stage rates in previous time step. + % T: Stage times in previous time step. + end + + function [a,b,c,s] = getTableau(obj) + [a,b,c,s] = obj.tsImplementation.getTableau(); + end + + % Returns quadrature weights for stages in one time step + function quadWeights = getTimeStepQuadrature(obj) + [~, b] = obj.getTableau(); + quadWeights = obj.k*b; + end + + % Use RK for first order form to step + function obj = step(obj) + obj.tsImplementation.step(); + [v, t] = obj.tsImplementation.getV(); + obj.w = v; + obj.t = t; + obj.n = obj.n + 1; + end + end + + methods (Static) + function k = getTimeStep(lambda, order) + default_arg('order', 4); + k = obj.tsImplementation.getTimeStep(lambda, order); + end + end + +end \ No newline at end of file
--- a/+time/SBPInTimeSecondOrderFormImplicit.m Sun Nov 04 12:36:30 2018 -0800 +++ b/+time/SBPInTimeSecondOrderFormImplicit.m Wed Nov 21 18:29:29 2018 -0800 @@ -19,11 +19,11 @@ default_arg('TYPE', []); default_arg('order', []); default_arg('blockSize',[]); - default_arg('do_scaling', true); + default_arg('do_scaling', false); m = length(v0); - default_arg('A', sparse(m, m)); + default_arg('A', speye(m, m)); default_arg('B', sparse(m, m)); default_arg('C', sparse(m, m));
--- a/.hgtags Sun Nov 04 12:36:30 2018 -0800 +++ b/.hgtags Wed Nov 21 18:29:29 2018 -0800 @@ -1,1 +1,4 @@ 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1 +0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1 +08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0 +08f3ffe63f484d02abce8df4df61e826f568193f Heimisson2018
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertIsMember.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,3 @@ +function assertIsMember(v, allowed) + assert(ismember(v, allowed), 'Expected ''%s'' to be in the set %s', inputname(1), toString(allowed)); +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertSize.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,16 @@ +% Assert that array A has the size s. +function assertSize(A,varargin) + if length(varargin) == 1 + s = varargin{1}; + errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), toString(s), toString(size(A))); + assert(all(size(A) == s), errmsg); + elseif length(varargin) == 2 + dim = varargin{1}; + s = varargin{2}; + + errmsg = sprintf('Expected %s to have size %d along dimension %d, got: %d',inputname(1), s, dim, size(A,dim)); + assert(size(A,dim) == s, errmsg); + else + error('Expected 2 or 3 arguments to assertSize()'); + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/assertStructFields.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,12 @@ +% Assert that the struct s has the all the field names in the cell array fns. +function assertStructFields(s, fns) + assertType(s, 'struct'); + assertType(fns, 'cell'); + + ok = ismember(fns, fieldnames(s)); + if ~all(ok) + str1 = sprintf("'%s' must have the fields %s\n", inputname(1), toString(fns)); + str2 = sprintf("The following fields are missing: %s", toString(fns(~ok))); + error(str1 + str2); + end +end
--- a/assert_size.m Sun Nov 04 12:36:30 2018 -0800 +++ b/assert_size.m Wed Nov 21 18:29:29 2018 -0800 @@ -1,16 +1,5 @@ % Assert that array A has the size s. function assert_size(A,s) - errmsg = sprintf('Expected %s to have size %s, got: %s',inputname(1), format_vector(s), format_vector(size(A))); - assert(all(size(A) == s),errmsg); -end - -function str = format_vector(a) - l = length(a); - str = sprintf('[%d',a(1)); - - for i = 2:l - str = [str sprintf(', %d',a(i))]; - end - - str = [str ']']; + warning('Use assertSize() instead!') + assertSize(A,s); end \ No newline at end of file
--- a/diracDiscr.m Sun Nov 04 12:36:30 2018 -0800 +++ b/diracDiscr.m Wed Nov 21 18:29:29 2018 -0800 @@ -1,4 +1,37 @@ -function ret = diracDiscr(x_0in , x , m_order, s_order, H, h) + +function d = diracDiscr(x_s, x, m_order, s_order, H) + % n-dimensional delta function + % x_s: source point coordinate vector, e.g. [x, y] or [x, y, z]. + % x: cell array of grid point column vectors for each dimension. + % m_order: Number of moment conditions + % s_order: Number of smoothness conditions + % H: cell array of 1D norm matrices + + dim = length(x_s); + d_1D = cell(dim,1); + + % If 1D, non-cell input is accepted + if dim == 1 && ~iscell(x) + d = diracDiscr1D(x_s, x, m_order, s_order, H); + + else + for i = 1:dim + d_1D{i} = diracDiscr1D(x_s(i), x{i}, m_order, s_order, H{i}); + end + + d = d_1D{dim}; + for i = dim-1: -1: 1 + % Perform outer product, transpose, and then turn into column vector + d = (d_1D{i}*d')'; + d = d(:); + end + end + +end + + +% Helper function for 1D delta functions +function ret = diracDiscr1D(x_0in , x , m_order, s_order, H) m = length(x); @@ -14,7 +47,7 @@ tot = m_order+s_order; S = []; M = []; - + % Get interior grid spacing middle = floor(m/2); h = x(middle+1) - x(middle); @@ -53,9 +86,9 @@ x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); norm = fnorm(poss)/h; index = poss; - + % Interior - else + else pol = (x(poss)-x(poss(1)))/(x(poss(end))-x(poss(1))); x_0 = (x_0in-x(poss(1)))/(x(poss(end))-x(poss(1))); norm = fnorm(poss)/h;
--- a/diracDiscrTest.m Sun Nov 04 12:36:30 2018 -0800 +++ b/diracDiscrTest.m Wed Nov 21 18:29:29 2018 -0800 @@ -10,7 +10,7 @@ for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test left boundary grid points x0s = xl + [0, h, 2*h]; @@ -33,11 +33,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test random points near left boundary x0s = xl + 2*h*rand(1,10); @@ -60,11 +60,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test right boundary grid points x0s = xr-[0, h, 2*h]; @@ -87,11 +87,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test random points near right boundary x0s = xr - 2*h*rand(1,10); @@ -114,11 +114,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test interior grid points m_half = round(m/2); @@ -142,11 +142,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test random points in interior x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); @@ -170,11 +170,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test points outisde grid x0s = [xl-1.1*h, xr+1.1*h]; @@ -197,11 +197,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test all grid points x0s = x; @@ -224,11 +224,11 @@ orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond); + [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond); % Test halfway between all grid points x0s = 1/2*( x(2:end)+x(1:end-1) ); @@ -247,81 +247,224 @@ end end -function testAllGPStaggered(testCase) +% function testAllGPStaggered(testCase) + +% orders = [2, 4, 6]; +% mom_conds = orders; + +% for o = 1:length(orders) +% order = orders(o); +% mom_cond = mom_conds(o); +% [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + +% % Test all grid points +% x0s = x; + +% for j = 1:length(fs) +% f = fs{j}; +% fx = f(x); +% for i = 1:length(x0s) +% x0 = x0s(i); +% delta = diracDiscr(x0, x, mom_cond, 0, H); +% integral = delta'*H*fx; +% err = abs(integral - f(x0)); +% testCase.verifyLessThan(err, 1e-12); +% end +% end +% end +% end + +% function testHalfGPStaggered(testCase) + +% orders = [2, 4, 6]; +% mom_conds = orders; + +% for o = 1:length(orders) +% order = orders(o); +% mom_cond = mom_conds(o); +% [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + +% % Test halfway between all grid points +% x0s = 1/2*( x(2:end)+x(1:end-1) ); + +% for j = 1:length(fs) +% f = fs{j}; +% fx = f(x); +% for i = 1:length(x0s) +% x0 = x0s(i); +% delta = diracDiscr(x0, x, mom_cond, 0, H); +% integral = delta'*H*fx; +% err = abs(integral - f(x0)); +% testCase.verifyLessThan(err, 1e-12); +% end +% end +% end +% end + +% function testRandomStaggered(testCase) + +% orders = [2, 4, 6]; +% mom_conds = orders; + +% for o = 1:length(orders) +% order = orders(o); +% mom_cond = mom_conds(o); +% [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + +% % Test random points within grid boundaries +% x0s = xl + (xr-xl)*rand(1,300); + +% for j = 1:length(fs) +% f = fs{j}; +% fx = f(x); +% for i = 1:length(x0s) +% x0 = x0s(i); +% delta = diracDiscr(x0, x, mom_cond, 0, H); +% integral = delta'*H*fx; +% err = abs(integral - f(x0)); +% testCase.verifyLessThan(err, 1e-12); +% end +% end +% end +% end + +%=============== 2D tests ============================== +function testAllGP2D(testCase) orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); % Test all grid points - x0s = x; + x0s = X; for j = 1:length(fs) f = fs{j}; - fx = f(x); + fx = f(X(:,1), X(:,2)); for i = 1:length(x0s) - x0 = x0s(i); + x0 = x0s(i,:); delta = diracDiscr(x0, x, mom_cond, 0, H); - integral = delta'*H*fx; - err = abs(integral - f(x0)); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2))); testCase.verifyLessThan(err, 1e-12); end end end end -function testHalfGPStaggered(testCase) +function testAllRandom2D(testCase) orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); - % Test halfway between all grid points - x0s = 1/2*( x(2:end)+x(1:end-1) ); + xl = xlims{1}; + xr = xlims{2}; + yl = ylims{1}; + yr = ylims{2}; + + % Test random points, even outside grid + Npoints = 100; + x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ... + (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1) ]; for j = 1:length(fs) f = fs{j}; - fx = f(x); + fx = f(X(:,1), X(:,2)); for i = 1:length(x0s) - x0 = x0s(i); + x0 = x0s(i,:); delta = diracDiscr(x0, x, mom_cond, 0, H); - integral = delta'*H*fx; - err = abs(integral - f(x0)); + integral = delta'*H_global*fx; + + % Integral should be 0 if point is outside grid + if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2))); + end testCase.verifyLessThan(err, 1e-12); end end end end -function testRandomStaggered(testCase) +%=============== 3D tests ============================== +function testAllGP3D(testCase) orders = [2, 4, 6]; mom_conds = orders; - + for o = 1:length(orders) order = orders(o); mom_cond = mom_conds(o); - [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond); + [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); - % Test random points within grid boundaries - x0s = xl + (xr-xl)*rand(1,300); + % Test all grid points + x0s = X; for j = 1:length(fs) f = fs{j}; - fx = f(x); + fx = f(X(:,1), X(:,2), X(:,3)); for i = 1:length(x0s) - x0 = x0s(i); + x0 = x0s(i,:); delta = diracDiscr(x0, x, mom_cond, 0, H); - integral = delta'*H*fx; - err = abs(integral - f(x0)); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2), x0(3))); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllRandom3D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + + xl = xlims{1}; + xr = xlims{2}; + yl = ylims{1}; + yr = ylims{2}; + zl = zlims{1}; + zr = zlims{2}; + + % Test random points, even outside grid + Npoints = 200; + x0s = [(xl-3*h{1}) + (xr-xl+6*h{1})*rand(Npoints,1), ... + (yl-3*h{2}) + (yr-yl+6*h{2})*rand(Npoints,1), ... + (zl-3*h{3}) + (zr-zl+6*h{3})*rand(Npoints,1) ]; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2), X(:,3)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H_global*fx; + + % Integral should be 0 if point is outside grid + if x0(1) < xl || x0(1) > xr || x0(2) < yl || x0(2) > yr || x0(3) < zl || x0(3) > zr + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2), x0(3))); + end testCase.verifyLessThan(err, 1e-12); end end @@ -329,8 +472,10 @@ end +% ====================================================== % ============== Setup functions ======================= -function [xl, xr, m, h, x, H, fs] = setupStuff(order, mom_cond) +% ====================================================== +function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond) % Grid xl = -3; @@ -353,6 +498,79 @@ end +function [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; + + m = [15, 16]; + g = grid.equidistant(m, xlims, ylims); + X = g.points(); + x = g.x; + + % Quadrature + opsx = sbp.D2Standard(m(1), xlims, order); + opsy = sbp.D2Standard(m(2), ylims, order); + Hx = opsx.H; + Hy = opsy.H; + H = {Hx, Hy}; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x,y) (x/Lx + y/Ly).^p; + end + + % Grid spacing in interior + mm = round(m/2); + hx = x{1}(mm(1)+1) - x{1}(mm(1)); + hy = x{2}(mm(2)+1) - x{2}(mm(2)); + h = {hx, hy}; + +end + +function [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond) + + % Grid + xlims = {-3, 20}; + ylims = {-11,5}; + zlims = {2,4}; + Lx = xlims{2} - xlims{1}; + Ly = ylims{2} - ylims{1}; + Lz = zlims{2} - zlims{1}; + + m = [13, 14, 15]; + g = grid.equidistant(m, xlims, ylims, zlims); + X = g.points(); + x = g.x; + + % Quadrature + opsx = sbp.D2Standard(m(1), xlims, order); + opsy = sbp.D2Standard(m(2), ylims, order); + opsz = sbp.D2Standard(m(3), zlims, order); + Hx = opsx.H; + Hy = opsy.H; + Hz = opsz.H; + H = {Hx, Hy, Hz}; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x,y,z) (x/Lx + y/Ly + z/Lz).^p; + end + + % Grid spacing in interior + mm = round(m/2); + hx = x{1}(mm(1)+1) - x{1}(mm(1)); + hy = x{2}(mm(2)+1) - x{2}(mm(2)); + hz = x{3}(mm(3)+1) - x{3}(mm(3)); + h = {hx, hy, hz}; + +end + function [xl, xr, m, h, x, H, fs] = setupStaggered(order, mom_cond) % Grid @@ -375,5 +593,3 @@ end end - -
--- a/mononomial.m Sun Nov 04 12:36:30 2018 -0800 +++ b/mononomial.m Wed Nov 21 18:29:29 2018 -0800 @@ -1,8 +1,17 @@ -function y = mononomial(x, k) - if k < 0 - y = x*0; +% calculate a N-D mononomial with powers k in points x: +% z = x(:,1).^k(1) * x(:,2).^k(2) * ... +function z = mononomial(x, k) + assert(size(x,2) == length(k), 'k must have the same length as the width of x'); + + if any(k < 0) + z = x(:,1)*0; return end - y = x.^k/factorial(k); + + denom = prod(factorial(k)); + + for i = 1:length(k) + x(:,i) = x(:,i).^k(i); + end + z = prod(x,2)/denom; end -
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/nextColor.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,5 @@ +function c = nextColor(ah) + default_arg('ah', gca); + + c = ah.ColorOrder(ah.ColorOrderIndex, :); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/pointIndex.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,4 @@ +% Get the index of the points p within the tall array of points ps +function [I, ok] = pointIndex(p, ps) + [ok, I] = ismember(p, ps, 'rows'); +end
--- a/spdiag.m Sun Nov 04 12:36:30 2018 -0800 +++ b/spdiag.m Wed Nov 21 18:29:29 2018 -0800 @@ -5,6 +5,6 @@ a = a'; end - n = length(a)-abs(i); + n = length(a)+abs(i); A = spdiags(a,i,n,n); end \ No newline at end of file
--- a/spdiagVariable.m Sun Nov 04 12:36:30 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,17 +0,0 @@ -function A = spdiagVariable(a,i) - default_arg('i',0); - - if isrow(a) - a = a'; - end - - n = length(a)+abs(i); - - if i > 0 - a = [sparse(i,1); a]; - elseif i < 0 - a = [a; sparse(abs(i),1)]; - end - - A = spdiags(a,i,n,n); -end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagsPeriodic.m Wed Nov 21 18:29:29 2018 -0800 @@ -0,0 +1,60 @@ +function A = spdiagsPeriodic(vals,diags) + % Creates an m x m periodic discretization matrix. + % vals - m x ndiags matrix of values + % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on + % vals that are not on main diagonal are going to spill over to + % off-diagonal corners. + + default_arg('diags',0); + + [m, ~] = size(vals); + + A = sparse(m,m); + + for i = 1:length(diags) + + d = diags(i); + a = vals(:,i); + + % Sub-diagonals + if d < 0 + a_bulk = a(1+abs(d):end); + a_corner = a(1:1+abs(d)-1); + corner_diag = m-abs(d); + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Super-diagonals + elseif d > 0 + a_bulk = a(1:end-d); + a_corner = a(end-d+1:end); + corner_diag = -m + d; + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Main diagonal + else + A = A + spdiagVariable(a, 0); + end + + end + +end + +function A = spdiagVariable(a,i) + default_arg('i',0); + + if isrow(a) + a = a'; + end + + n = length(a)+abs(i); + + if i > 0 + a = [sparse(i,1); a]; + elseif i < 0 + a = [a; sparse(abs(i),1)]; + end + + A = spdiags(a,i,n,n); +end
--- a/spdiagsVariablePeriodic.m Sun Nov 04 12:36:30 2018 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,42 +0,0 @@ -function A = spdiagsVariablePeriodic(vals,diags) - % Creates an m x m periodic discretization matrix. - % vals - m x ndiags matrix of values - % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on - % vals that are not on main diagonal are going to spill over to - % off-diagonal corners. - - default_arg('diags',0); - - [m, ~] = size(vals); - - A = sparse(m,m); - - for i = 1:length(diags) - - d = diags(i); - a = vals(:,i); - - % Sub-diagonals - if d < 0 - a_bulk = a(1+abs(d):end); - a_corner = a(1:1+abs(d)-1); - corner_diag = m-abs(d); - A = A + spdiagVariable(a_bulk, d); - A = A + spdiagVariable(a_corner, corner_diag); - - % Super-diagonals - elseif d > 0 - a_bulk = a(1:end-d); - a_corner = a(end-d+1:end); - corner_diag = -m + d; - A = A + spdiagVariable(a_bulk, d); - A = A + spdiagVariable(a_corner, corner_diag); - - % Main diagonal - else - A = A + spdiagVariable(a, 0); - end - - end - -end \ No newline at end of file
--- a/stripeMatrixPeriodic.m Sun Nov 04 12:36:30 2018 -0800 +++ b/stripeMatrixPeriodic.m Wed Nov 21 18:29:29 2018 -0800 @@ -1,8 +1,8 @@ -% Creates a periodic discretization matrix of size n x n +% Creates a periodic discretization matrix of size n x n % with the values of val on the diagonals diag. % A = stripeMatrix(val,diags,n) function A = stripeMatrixPeriodic(val,diags,n) D = ones(n,1)*val; - A = spdiagsVariablePeriodic(D,diags); + A = spdiagsPeriodic(D,diags); end \ No newline at end of file
--- a/vandermonde.m Sun Nov 04 12:36:30 2018 -0800 +++ b/vandermonde.m Wed Nov 21 18:29:29 2018 -0800 @@ -1,10 +1,15 @@ % Create vandermonde matrix for points x and polynomials of order p -% x and p are vectors -% v is a length(x) by length(p) matrix +% x is a list of N points of size [N,dim], +% p is a list of polynomial orders of size [M, dim]. +% the given mononomials are evaluated and the NxM matrix V is returned. function V = vandermonde(x, p) - V = sym(zeros(length(x), length(p))); % Is there a way to make this work for both double and sym + assert(size(x,2) == size(p,2), 'x and p must have the same number of columns') + n = size(x,1); + m = size(p,1); - for i = 1:length(p) - V(:, i) = mononomial(x,p(i)); + for i = 1:m + V(:,i) = mononomial(x, p(i,:)); end + + assertSize(V,[n,m]); end