Mercurial > repos > public > sbplib
changeset 1071:92cb03e64ca4 feature/grids/LaplaceSquared
Merge with default
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Thu, 20 Sep 2018 12:05:20 +0200 |
| parents | f6b3af6febf3 (current diff) 006defd0247b (diff) |
| children | 6468a5f6ec79 |
| files | +multiblock/multiblockgrid.m +multiblock/stitchSchemes.m |
| diffstat | 26 files changed, 2298 insertions(+), 204 deletions(-) [+] |
line wrap: on
line diff
--- a/+grid/Cartesian.m Tue Jul 24 21:08:08 2018 -0700 +++ b/+grid/Cartesian.m Thu Sep 20 12:05:20 2018 +0200 @@ -5,6 +5,7 @@ m % Number of points in each direction x % Cell array of vectors with node placement for each dimension. h % Spacing/Scaling + lim % Cell array of left and right boundaries for each dimension. end % General d dimensional grid with n points @@ -27,6 +28,11 @@ end obj.h = []; + + 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 Tue Jul 24 21:08:08 2018 -0700 +++ b/+grid/evalOn.m Thu Sep 20 12:05:20 2018 +0200 @@ -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/DiffOp.m Tue Jul 24 21:08:08 2018 -0700 +++ b/+multiblock/DiffOp.m Thu Sep 20 12:05:20 2018 +0200 @@ -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
--- a/+multiblock/Grid.m Tue Jul 24 21:08:08 2018 -0700 +++ b/+multiblock/Grid.m Thu Sep 20 12:05:20 2018 +0200 @@ -77,7 +77,7 @@ % Collect number of points in each block N = zeros(1,nBlocks); for i = 1:nBlocks - N(i) = obj.grids{i}.N(); + N(i) = obj.grids{i}.N()*nComponents; end gfs = blockmatrix.fromMatrix(gf, {N,1});
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+multiblock/evalOn.m Thu Sep 20 12:05:20 2018 +0200 @@ -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 Tue Jul 24 21:08:08 2018 -0700 +++ /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 Tue Jul 24 21:08:08 2018 -0700 +++ /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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_2.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,50 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<3) + error(['Operator requires at least ' num2str(3) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + % D1 operator + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + scheme_width = 3; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + Mm1 = -c(r-1)/2 - c(r)/2; + M0 = c(r-1)/2 + c(r) + c(r+1)/2; + Mp1 = -c(r)/2 - c(r+1)/2; + + vals = [Mm1,M0,Mp1]; + diags = -scheme_radius : scheme_radius; + M = spdiagsPeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_4.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,57 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_4(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<5) + error(['Operator requires at least ' num2str(5) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + S = d1_l*d1_l' + d1_r*d1_r'; + + % D1 operator + stencil = [1/12 -2/3 0 2/3 -1/12]; + diags = -2:2; + Q = stripeMatrixPeriodic(stencil, diags, m); + D1 = HI*(Q - 1/2*e_l*e_l' + 1/2*e_r*e_r'); + + + scheme_width = 5; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + % Note: these coefficients are for -M. + Mm2 = -1/8*c(r-2) + 1/6*c(r-1) - 1/8*c(r); + Mm1 = 1/6 *c(r-2) + 1/2*c(r-1) + 1/2*c(r) + 1/6*c(r+1); + M0 = -1/24*c(r-2)- 5/6*c(r-1) - 3/4*c(r) - 5/6*c(r+1) - 1/24*c(r+2); + Mp1 = 0 * c(r-2) + 1/6*c(r-1) + 1/2*c(r) + 1/2*c(r+1) + 1/6 *c(r+2); + Mp2 = 0 * c(r-2) + 0 * c(r-1) - 1/8*c(r) + 1/6*c(r+1) - 1/8 *c(r+2); + + vals = -[Mm2,Mm1,M0,Mp1,Mp2]; + diags = -scheme_radius : scheme_radius; + M = spdiagsPeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + + end + D2 = @D2_fun; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_6.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,58 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_6(m,h) + % m = number of unique grid points, i.e. h = L/m; + + if(m<7) + error(['Operator requires at least ' num2str(7) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + + % D1 operator + diags = -3:3; + stencil = [-1/60 9/60 -45/60 0 45/60 -9/60 1/60]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + % D2 operator + scheme_width = 7; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + 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); + Mm1 = -c(r-3)/0.40e2 - 0.3e1/0.10e2 * c(r-2) - 0.3e1/0.10e2 * c(r+1) - c(r+2)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r-1); + M0 = c(r-3)/0.180e3 + c(r-2)/0.8e1 + 0.19e2/0.20e2 * c(r-1) + 0.19e2/0.20e2 * c(r+1) + c(r+2)/0.8e1 + c(r+3)/0.180e3 + 0.101e3/0.180e3 * c(r); + Mp1 = -c(r-2)/0.40e2 - 0.3e1/0.10e2 * c(r-1) - 0.3e1/0.10e2 * c(r+2) - c(r+3)/0.40e2 - 0.17e2/0.40e2 * c(r) - 0.17e2/0.40e2 * c(r+1); + Mp2 = c(r-1)/0.20e2 - 0.3e1/0.10e2 * c(r+1) + c(r+3)/0.20e2 + 0.7e1/0.40e2 * c(r) + 0.7e1/0.40e2 * c(r+2); + Mp3 = c(r+1)/0.40e2 + c(r+2)/0.40e2 - 0.11e2/0.360e3 * c(r) - 0.11e2/0.360e3 * c(r+3); + + vals = [Mm3,Mm2,Mm1,M0,Mp1,Mp2,Mp3]; + diags = -scheme_radius : scheme_radius; + M = spdiagsPeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; + + +end
--- a/+sbp/D2Variable.m Tue Jul 24 21:08:08 2018 -0700 +++ b/+sbp/D2Variable.m Thu Sep 20 12:05:20 2018 +0200 @@ -26,22 +26,39 @@ obj.x = linspace(x_l,x_r,m)'; switch order + + case 6 + + [obj.H, obj.HI, obj.D1, obj.D2, ... + ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... + obj.d1_l, obj.d1_r] = ... + sbp.implementations.d4_variable_6(m, obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + case 4 [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_4(m,obj.h); obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R case 2 [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... obj.e_r, obj.d1_l, obj.d1_r] = ... sbp.implementations.d2_variable_2(m,obj.h); obj.borrowing.M.d1 = 0.3636363636; % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R otherwise error('Invalid operator order %d.',order); end - + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, obj.m = m; obj.M = []; end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2VariablePeriodic.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,71 @@ +classdef D2VariablePeriodic < sbp.OpSet + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + D2 % SBP operator for second derivative + M % Norm matrix, second derivative + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D2VariablePeriodic(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + obj.h = L/m; + x = linspace(x_l,x_r,m+1)'; + obj.x = x(1:end-1); + + switch order + + case 6 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_6(m,obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + + case 4 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_4(m,obj.h); + obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R + case 2 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_2(m,obj.h); + obj.borrowing.M.d1 = 0.3636363636; + % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R + + otherwise + error('Invalid operator order %d.',order); + end + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, + + obj.m = m; + obj.M = []; + end + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end + + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Elastic2dCurvilinear.m Thu Sep 20 12:05:20 2018 +0200 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Elastic2dVariable.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,512 @@ +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 +% 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 + + 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 + + % Borrowing from H, M, and R + thH + thM + thR + + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + E % E{i}^T picks out component i + + H_boundary % Boundary inner products + + % Kroneckered norms and coefficients + RHOi_kron + Hi_kron + end + + methods + + function obj = Elastic2dVariable(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; + + assert(isa(g, 'grid.Cartesian')) + + lambda = grid.evalOn(g, lambda_fun); + mu = grid.evalOn(g, mu_fun); + rho = grid.evalOn(g, rho_fun); + 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); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, 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; + + % Better names + obj.thR{i} = ops{i}.borrowing.R.delta_D; + obj.thM{i} = ops{i}.borrowing.M.d1; + obj.thH{i} = ops{i}.borrowing.H11; + 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 ======== + LAMBDA = spdiag(lambda); + obj.LAMBDA = LAMBDA; + MU = spdiag(mu); + obj.MU = MU; + RHO = spdiag(rho); + obj.RHO = RHO; + obj.RHOi = inv(RHO); + + obj.D1 = cell(dim,1); + obj.D2_lambda = cell(dim,1); + obj.D2_mu = 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}); + + % 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 + obj.D2_lambda{i} = sparse(m_tot); + obj.D2_mu{i} = sparse(m_tot); + end + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m(2) + D_lambda = D2{1}(lambda(ind(:,i))); + D_mu = D2{1}(mu(ind(:,i))); + + p = ind(:,i); + obj.D2_lambda{1}(p,p) = D_lambda; + obj.D2_mu{1}(p,p) = D_mu; + end + + for i = 1:m(1) + D_lambda = D2{2}(lambda(ind(i,:))); + D_mu = D2{2}(mu(ind(i,:))); + + p = ind(i,:); + obj.D2_lambda{2}(p,p) = D_lambda; + obj.D2_mu{2}(p,p) = D_mu; + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary = cell(dim,1); + obj.H_boundary{1} = H{2}; + obj.H_boundary{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 + D = D + E{i}*inv(RHO)*( d(i,j)*D2_lambda{i}*E{j}' +... + db(i,j)*D1{i}*LAMBDA*D1{j}*E{j}' ... + ); + D = D + E{i}*inv(RHO)*( d(i,j)*D2_mu{i}*E{j}' +... + db(i,j)*D1{j}*MU*D1{i}*E{j}' + ... + D2_mu{j}*E{i}' ... + ); + 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; + D1 = obj.D1; + + % 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); + for k = 1:dim + T_l{j}{i,k} = ... + -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... + -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... + -d(i,k)*MU*e_l{j}*d1_l{j}'; + + T_r{j}{i,k} = ... + d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... + +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... + +d(i,k)*MU*e_r{j}*d1_r{j}'; + + 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.Hi_kron = kron(obj.Hi, I_dim); + + % Misc. + obj.m = m; + obj.h = h; + 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; + + 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)); + + 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*B'*e*H_gamma*(e'*E{k}' ); + penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; + end + + % Free boundary condition + case {'F','f','Free','free','traction','Traction','t','T'} + closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); + penalty = penalty + E{k}*RHOi*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. + 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 + h_u = obj.h(j); + thR_u = obj.thR{j}*h_u; + thM_u = obj.thM{j}*h_u; + thH_u = obj.thH{j}*h_u; + + h_v = neighbour_scheme.h(j_v); + thR_v = neighbour_scheme.thR{j_v}*h_v; + thH_v = neighbour_scheme.thH{j_v}*h_v; + thM_v = neighbour_scheme.thM{j_v}*h_v; + + % alpha = penalty strength for normal component, beta for tangential + alpha_u = dim*lambda_u/(4*thH_u) + lambda_u/(4*thR_u) + mu_u/(2*thM_u); + alpha_v = dim*lambda_v/(4*thH_v) + lambda_v/(4*thR_v) + mu_v/(2*thM_v); + beta_u = mu_u/(2*thH_u) + mu_u/(4*thR_u); + beta_v = mu_v/(2*thH_v) + mu_v/(4*thR_v); + alpha = alpha_u + alpha_v; + beta = beta_u + beta_v; + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + strength = @(i,j) tuning*(d(i,j)*alpha + db(i,j)*beta); + + % Preallocate + closure = sparse(dim*m_tot_u, dim*m_tot_u); + penalty = sparse(dim*m_tot_u, dim*m_tot_v); + + % Loop over components that penalties end up on + for i = 1:dim + closure = closure - E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e'*E{i}'; + penalty = penalty + E{i}*RHOi*Hi*e*strength(i,j)*H_gamma*e_v'*E_v{i}'; + + closure = closure - 1/2*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' + 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 + otherwise + error(['No such operator: operator = ' op{i}]); + end + end + + end + + function N = size(obj) + N = obj.dim*prod(obj.m); + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Heat2dCurvilinear.m Thu Sep 20 12:05:20 2018 +0200 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Heat2dVariable.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,276 @@ +classdef Heat2dVariable < scheme.Scheme + +% Discretizes the Laplacian with variable coefficent, +% 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 + + H_boundary % Boundary inner products + + end + + methods + + function obj = Heat2dVariable(g ,order, kappa_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('kappa_fun', @(x,y) 0*x+1); + dim = 2; + + assert(isa(g, 'grid.Cartesian')) + + kappa = grid.evalOn(g, kappa_fun); + m = g.size(); + m_tot = g.N(); + + h = g.scaling(); + lim = g.lim; + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, 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; + + 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}); + + % 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 + obj.D2_kappa{i} = sparse(m_tot); + end + ind = grid.funcToMatrix(g, 1:m_tot); + + for i = 1:m(2) + D_kappa = D2{1}(kappa(ind(:,i))); + p = ind(:,i); + obj.D2_kappa{1}(p,p) = D_kappa; + end + + for i = 1:m(1) + D_kappa = D2{2}(kappa(ind(i,:))); + p = ind(i,:); + obj.D2_kappa{2}(p,p) = D_kappa; + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary = cell(dim,1); + obj.H_boundary{1} = H{2}; + obj.H_boundary{2} = H{1}; + + % Differentiation matrix D (without SAT) + D2_kappa = obj.D2_kappa; + D1 = obj.D1; + D = sparse(m_tot,m_tot); + for i = 1:dim + D = D + D2_kappa{i}; + end + obj.D = D; + %=========================================% + + % Misc. + obj.m = m; + obj.h = h; + 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; + d = obj.d1_r; + case -1 + e = obj.e_l; + d = obj.d1_l; + end + + Hi = obj.Hi; + 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 = Hi*e{j}*kappa_gamma*H_gamma; + % penalty is for normal derivative and not for derivative, hence the sign. + + % 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/bcSetup.m Tue Jul 24 21:08:08 2018 -0700 +++ b/+scheme/bcSetup.m Thu Sep 20 12:05:20 2018 +0200 @@ -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]; - assertSize(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/.hgtags Tue Jul 24 21:08:08 2018 -0700 +++ b/.hgtags Thu Sep 20 12:05:20 2018 +0200 @@ -1,2 +1,4 @@ 18c023aaf3f79cbe2b9b1cf547d80babdaa1637d v0.1 0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1 +b723495cdb2f96314d7b3f0aa79723a7dc088c7d v0.2 +08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/LICENSE.txt Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,25 @@ +MIT License + +Copyright (c) +2015-2018 Jonatan Werpers +2015-2018 Martin Almquist +2016-2018 Ylva Rydin +2018 Vidar Stiernström + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/README.md Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,2 @@ +# SBPLIB +sbplib is a library of primitives and help functions for working with summation-by-parts finite differences in Matlab. To use sbplib download the code and add the sbplib folder to the matlab path.
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/kroneckerDelta.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,6 @@ +function d = kroneckerDelta(i,j) + +d = 0; +if i==j + d = 1; +end \ No newline at end of file
--- a/spdiag.m Tue Jul 24 21:08:08 2018 -0700 +++ b/spdiag.m Thu Sep 20 12:05:20 2018 +0200 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagsPeriodic.m Thu Sep 20 12:05:20 2018 +0200 @@ -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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stripeMatrixPeriodic.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,8 @@ +% 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 = spdiagsPeriodic(D,diags); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/structArray.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,16 @@ +% % Usage example: +% c = structArray({'a','b'}, { +% 1, 2; +% 3, 4; +% }); + +function c = structArray(fields, values) + assert(length(fields) == size(values, 2), 'Number of fields and number of colums of ''values'' must be equal'); + c = struct(); + + for i = 1:size(values, 1) + for j = 1:length(fields) + c(i).(fields{j}) = values{i,j}; + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/structCellArray.m Thu Sep 20 12:05:20 2018 +0200 @@ -0,0 +1,16 @@ +% % Usage example: +% c = structCellArray({'a','b'}, { +% 1, 2; +% 3, 4; +% }); + +function c = structCellArray(fields, values) + assert(length(fields) == size(values, 2), 'Number of fields and number of colums of ''values'' must be equal'); + c = cell(1, size(values, 1)); + + for i = 1:size(values, 1) + for j = 1:length(fields) + c{i}.(fields{j}) = values{i,j}; + end + end +end