Mercurial > repos > public > sbplib
changeset 968:a4ad90b37998 feature/poroelastic
Merge with default.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sun, 23 Dec 2018 14:39:31 +0100 |
| parents | 368a2773f78b (diff) 2412f407749a (current diff) |
| children | adae8063ea2f |
| files | +multiblock/DiffOp.m +sbp/+implementations/intOpAWW_orders_2to2_ratio2to1.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F1_accF2C2.m +sbp/+implementations/intOpAWW_orders_2to2_ratio_2to1_accC2F2_accF2C1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio2to1.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F2_accF2C3.m +sbp/+implementations/intOpAWW_orders_4to4_ratio_2to1_accC2F3_accF2C2.m +sbp/+implementations/intOpAWW_orders_6to6_ratio2to1.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F3_accF2C4.m +sbp/+implementations/intOpAWW_orders_6to6_ratio_2to1_accC2F4_accF2C3.m +sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5.m +sbp/+implementations/intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4.m +sbp/InterpAWW.m +sbp/InterpMC.m +scheme/Elastic2dVariable.m +scheme/Wave.m .hgtags |
| diffstat | 12 files changed, 1297 insertions(+), 94 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+multiblock/DiffOp.m Sun Dec 23 14:39:31 2018 +0100 @@ -150,6 +150,74 @@ end end + % Get a boundary operator specified by opName for the given boundary/BoundaryGroup + function op = getBoundaryOperatorWrapper(obj, opName, boundary, blockmatrixDiv) + default_arg('blockmatrixDiv', obj.blockmatrixDiv{1}); + + switch class(boundary) + case 'cell' + blockId = boundary{1}; + localOp = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2}); + + div = {blockmatrixDiv, size(localOp,2)}; + blockOp = blockmatrix.zero(div); + blockOp{blockId,1} = localOp; + op = blockmatrix.toMatrix(blockOp); + return + case 'multiblock.BoundaryGroup' + op = sparse(sum(blockmatrixDiv),0); + for i = 1:length(boundary) + op = [op, obj.getBoundaryOperatorWrapper(opName, boundary{i}, blockmatrixDiv)]; + end + otherwise + error('Unknown boundary indentifier') + end + end + + % Get a boundary cell of operators, specified by opName for the given boundary/BoundaryGroup + function opCell = getBoundaryCellOperator(obj, opName, boundary, blockmatrixDiv) + default_arg('blockmatrixDiv', obj.blockmatrixDiv{1}); + + % Get size of cell + switch class(boundary) + case 'cell' + blockId = boundary{1}; + localCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2}); + case 'multiblock.BoundaryGroup' + blockId = boundary{1}{1}; + localCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{1}{2}); + otherwise + error('Unknown boundary indentifier') + end + + % Loop over cell elements and build the boundary operator in each cell + opCell = cell(size(localCell)); + for i = 1:numel(opCell) + switch class(boundary) + case 'cell' + blockId = boundary{1}; + localOpCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2}); + localOp = localOpCell{i}; + + div = {blockmatrixDiv, size(localOp,2)}; + blockOp = blockmatrix.zero(div); + blockOp{blockId,1} = localOp; + op = blockmatrix.toMatrix(blockOp); + opCell{i} = op; + + case 'multiblock.BoundaryGroup' + op = sparse(sum(blockmatrixDiv),0); + for j = 1:length(boundary) + localCell = obj.getBoundaryCellOperator(opName, boundary{j}, blockmatrixDiv); + op = [op, localCell{i}]; + end + opCell{i} = op; + otherwise + error('Unknown boundary indentifier') + end + end + end + function op = getBoundaryQuadrature(obj, boundary) opName = 'H'; switch class(boundary)
--- a/+multiblock/Grid.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+multiblock/Grid.m Sun Dec 23 14:39:31 2018 +0100 @@ -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/+sbp/+implementations/d2_variable_2.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+sbp/+implementations/d2_variable_2.m Sun Dec 23 14:39:31 2018 +0100 @@ -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 Wed Dec 12 23:16:44 2018 +0100 +++ b/+sbp/+implementations/d2_variable_4.m Sun Dec 23 14:39:31 2018 +0100 @@ -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/d4_variable_6.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+sbp/+implementations/d4_variable_6.m Sun Dec 23 14:39:31 2018 +0100 @@ -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/+scheme/Elastic2dVariable.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+scheme/Elastic2dVariable.m Sun Dec 23 14:39:31 2018 +0100 @@ -30,42 +30,52 @@ 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 + e_w, e_e, e_s, e_n - % 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 + H_w, H_e, H_s, H_n % Kroneckered norms and coefficients RHOi_kron Hi_kron + + % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. + theta_R % Borrowing (d1- D1)^2 from R + theta_H % First entry in norm matrix + theta_M % Borrowing d1^2 from M. + + % Structures used for adjoint optimization + B end methods - function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet) + % The coefficients can either be function handles or grid functions + function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet) default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); - default_arg('lambda_fun', @(x,y) 0*x+1); - default_arg('mu_fun', @(x,y) 0*x+1); - default_arg('rho_fun', @(x,y) 0*x+1); + default_arg('lambda', @(x,y) 0*x+1); + default_arg('mu', @(x,y) 0*x+1); + default_arg('rho', @(x,y) 0*x+1); dim = 2; assert(isa(g, 'grid.Cartesian')) - lambda = grid.evalOn(g, lambda_fun); - mu = grid.evalOn(g, mu_fun); - rho = grid.evalOn(g, rho_fun); + if isa(lambda, 'function_handle') + lambda = grid.evalOn(g, lambda); + end + if isa(mu, 'function_handle') + mu = grid.evalOn(g, mu); + end + if isa(rho, 'function_handle') + rho = grid.evalOn(g, rho); + end + m = g.size(); m_tot = g.N(); @@ -87,15 +97,9 @@ % Borrowing constants for i = 1:dim - beta = ops{i}.borrowing.R.delta_D; - obj.H11{i} = ops{i}.borrowing.H11; - obj.phi{i} = beta/obj.H11{i}; - obj.gamma{i} = ops{i}.borrowing.M.d1; - - % Better names - obj.thR{i} = ops{i}.borrowing.R.delta_D; - obj.thM{i} = ops{i}.borrowing.M.d1; - obj.thH{i} = ops{i}.borrowing.H11; + obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D; + obj.theta_H{i} = h(i)*ops{i}.borrowing.H11; + obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1; end I = cell(dim,1); @@ -146,6 +150,10 @@ 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.e_w = obj.e_l{1}; + obj.e_e = obj.e_r{1}; + obj.e_s = obj.e_l{2}; + obj.e_n = obj.e_r{2}; obj.d1_l{1} = kron(d1_l{1},I{2}); obj.d1_l{2} = kron(I{1},d1_l{2}); @@ -183,6 +191,11 @@ 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}}; + obj.H_w = H{2}; + obj.H_e = H{2}; + obj.H_s = H{1}; + obj.H_n = H{1}; % E{i}^T picks out component i. E = cell(dim,1); @@ -213,7 +226,7 @@ end end obj.D = D; - %=========================================% + %=========================================%' % Numerical traction operators for BC. % Because d1 =/= e0^T*D1, the numerical tractions are different @@ -237,20 +250,28 @@ tau_l{j} = cell(dim,1); tau_r{j} = cell(dim,1); + LAMBDA_l = e_l{j}'*LAMBDA*e_l{j}; + LAMBDA_r = e_r{j}'*LAMBDA*e_r{j}; + MU_l = e_l{j}'*MU*e_l{j}; + MU_r = e_r{j}'*MU*e_r{j}; + + [~, n_l] = size(e_l{j}); + [~, n_r] = size(e_r{j}); + % Loop over components for i = 1:dim - tau_l{j}{i} = sparse(m_tot,dim*m_tot); - tau_r{j}{i} = sparse(m_tot,dim*m_tot); + tau_l{j}{i} = sparse(n_l, dim*m_tot); + tau_r{j}{i} = sparse(n_r, dim*m_tot); for k = 1:dim T_l{j}{i,k} = ... - -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})... - -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})... - -d(i,k)*MU*e_l{j}*d1_l{j}'; + -d(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})... + -d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})... + -d(i,k)*MU_l*d1_l{j}'; T_r{j}{i,k} = ... - d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})... - +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})... - +d(i,k)*MU*e_r{j}*d1_r{j}'; + d(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})... + +d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})... + +d(i,k)*MU_r*d1_r{j}'; tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}'; tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}'; @@ -258,6 +279,19 @@ end end + + % Transpose T and tau to match boundary operator convention + for i = 1:dim + for j = 1:dim + tau_l{i}{j} = transpose(tau_l{i}{j}); + tau_r{i}{j} = transpose(tau_r{i}{j}); + for k = 1:dim + T_l{i}{j,k} = transpose(T_l{i}{j,k}); + T_r{i}{j,k} = transpose(T_r{i}{j,k}); + end + end + end + obj.T_l = T_l; obj.T_r = T_r; obj.tau_l = tau_l; @@ -275,6 +309,44 @@ obj.grid = g; obj.dim = dim; + % B, used for adjoint optimization + B = cell(dim, 1); + for i = 1:dim + B{i} = cell(m_tot, 1); + end + + for i = 1:dim + for j = 1:m_tot + B{i}{j} = sparse(m_tot, m_tot); + end + end + + ind = grid.funcToMatrix(g, 1:m_tot); + + % Direction 1 + for k = 1:m(1) + c = sparse(m(1),1); + c(k) = 1; + [~, B_1D] = ops{1}.D2(c); + for l = 1:m(2) + p = ind(:,l); + B{1}{(k-1)*m(2) + l}(p, p) = B_1D; + end + end + + % Direction 2 + for k = 1:m(2) + c = sparse(m(2),1); + c(k) = 1; + [~, B_1D] = ops{2}.D2(c); + for l = 1:m(1) + p = ind(l,:); + B{2}{(l-1)*m(2) + k}(p, p) = B_1D; + end + end + + obj.B = B; + end @@ -297,6 +369,7 @@ 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; @@ -316,33 +389,20 @@ % Dirichlet boundary condition case {'D','d','dirichlet','Dirichlet'} - phi = obj.phi{j}; - h = obj.h(j); - h11 = obj.H11{j}*h; - gamma = obj.gamma{j}; - - a_lambda = dim/h11 + 1/(h11*phi); - a_mu_i = 2/(gamma*h); - a_mu_ij = 2/h11 + 1/(h11*phi); - - d = @kroneckerDelta; % Kronecker delta - db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta - alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... - + d(i,j)* a_mu_i*MU ... - + db(i,j)*a_mu_ij*MU ); + alpha = obj.get_boundary_operator('alpha', boundary); % Loop over components that Dirichlet penalties end up on for i = 1:dim - C = T{k,i}; - A = -d(i,k)*alpha(i,j); - B = A + C; + C = transpose(T{k,i}); + A = -e*transpose(alpha{i,k}); + B = A + e*C; closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' ); penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma; end % Free boundary condition case {'F','f','Free','free','traction','Traction','t','T'} - closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} ); + closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}'; penalty = penalty + E{k}*RHOi*Hi*e*H_gamma; % Unknown boundary condition @@ -351,11 +411,32 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % type Struct that specifies the interface coupling. + % Fields: + % -- tuning: penalty strength, defaults to 1.2 + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.tuning = 1.2; + defaultType.interpolation = 'none'; + default_struct('type', defaultType); + + switch type.interpolation + case {'none', ''} + [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); + case {'op','OP'} + [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); + otherwise + error('Unknown type of interpolation: %s ', type.interpolation); + end + end + + function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) + tuning = type.tuning; + % u denotes the solution in the own domain % v denotes the solution in the neighbour domain % Operators without subscripts are from the own domain. - tuning = 1.2; % j is the coordinate direction of the boundary j = obj.get_boundary_number(boundary); @@ -387,27 +468,27 @@ %------------------------- % 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; + 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}; - 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; + 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 = 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; + [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 = alpha_ii_u + alpha_ii_v; + sigma_ij = alpha_ij_u + alpha_ij_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); + 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); @@ -415,20 +496,24 @@ % 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 - 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}; + closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*tau{i}'; + penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*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}'; + closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e'*E{k}'; + penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}*H_gamma*e_v'*E_v{k}'; end end end + function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) + error('Non-conforming interfaces not implemented yet.'); + end + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. function [j, nj] = get_boundary_number(obj, boundary) @@ -466,40 +551,76 @@ op = {op}; end - for i = 1:length(op) - switch op{i} + for k = 1:length(op) + switch op{k} case 'e' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.e_l{j}; + varargout{k} = obj.e_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.e_r{j}; + varargout{k} = obj.e_r{j}; end case 'd' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.d1_l{j}; + varargout{k} = obj.d1_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.d1_r{j}; + varargout{k} = obj.d1_r{j}; end - case 'H' - varargout{i} = obj.H_boundary{j}; case 'T' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.T_l{j}; + varargout{k} = obj.T_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.T_r{j}; + varargout{k} = obj.T_r{j}; end case 'tau' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.tau_l{j}; + varargout{k} = obj.tau_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + varargout{k} = obj.tau_r{j}; + end + case 'H' + varargout{k} = obj.H_boundary{j}; + case 'alpha' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + switch boundary + case {'w','W','west','West','s','S','south','South'} + e = obj.e_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.tau_r{j}; + e = obj.e_r{j}; end + + tuning = 1.2; + LAMBDA = obj.LAMBDA; + MU = obj.MU; + + dim = obj.dim; + theta_R = obj.theta_R{j}; + theta_H = obj.theta_H{j}; + theta_M = obj.theta_M{j}; + + a_lambda = dim/theta_H + 1/theta_R; + a_mu_i = 2/theta_M; + a_mu_ij = 2/theta_H + 1/theta_R; + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + alpha = cell(obj.dim, obj.dim); + + alpha_func = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ... + + d(i,j)* a_mu_i*MU ... + + db(i,j)*a_mu_ij*MU ); + for i = 1:obj.dim + for l = 1:obj.dim + alpha{i,l} = d(i,l)*alpha_func(i,j)*e; + end + end + + varargout{k} = alpha; otherwise - error(['No such operator: operator = ' op{i}]); + error(['No such operator: operator = ' op{k}]); end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/+rkparameters/rk4.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,12 @@ +function [a,b,c,s] = rk4() + +% Butcher tableau for classical RK$ +s = 4; +a = sparse(s,s); +a(2,1) = 1/2; +a(3,2) = 1/2; +a(4,3) = 1; +b = 1/6*[1; 2; 2; 1]; +c = [0; 1/2; 1/2; 1]; + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/ExplicitRungeKuttaDiscreteData.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,121 @@ +classdef ExplicitRungeKuttaDiscreteData < time.Timestepper + properties + D + S % Function handle for time-dependent data + data % Matrix of data vectors, one column per stage + k + t + v + m + n + order + a, b, c, s % Butcher tableau + K % Stage rates + U % Stage approximations + T % Stage times + end + + + methods + function obj = ExplicitRungeKuttaDiscreteData(D, S, data, k, t0, v0, order) + default_arg('order', 4); + default_arg('S', []); + default_arg('data', []); + + obj.D = D; + obj.S = S; + obj.k = k; + obj.t = t0; + obj.v = v0; + obj.m = length(v0); + obj.n = 0; + obj.order = order; + obj.data = data; + + switch order + case 4 + [obj.a, obj.b, obj.c, obj.s] = time.rkparameters.rk4(); + otherwise + error('That RK method is not available'); + end + + obj.K = sparse(obj.m, obj.s); + obj.U = sparse(obj.m, obj.s); + + end + + function [v,t,U,T,K] = getV(obj) + v = obj.v; + t = obj.t; + U = obj.U; % Stage approximations in previous time step. + T = obj.T; % Stage times in previous time step. + K = obj.K; % Stage rates in previous time step. + end + + function [a,b,c,s] = getTableau(obj) + a = obj.a; + b = obj.b; + c = obj.c; + 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) + v = obj.v; + a = obj.a; + b = obj.b; + c = obj.c; + s = obj.s; + S = obj.S; + dt = obj.k; + K = obj.K; + U = obj.U; + D = obj.D; + data = obj.data; + + for i = 1:s + U(:,i) = v; + for j = 1:i-1 + U(:,i) = U(:,i) + dt*a(i,j)*K(:,j); + end + + K(:,i) = D*U(:,i); + obj.T(i) = obj.t + c(i)*dt; + + % Data from continuous function and discrete time-points. + if ~isempty(S) + K(:,i) = K(:,i) + S(obj.T(i)); + end + if ~isempty(data) + K(:,i) = K(:,i) + data(:,obj.n*s + i); + end + + end + + obj.v = v + dt*K*b; + obj.t = obj.t + dt; + obj.n = obj.n + 1; + obj.U = U; + obj.K = K; + end + end + + + methods (Static) + function k = getTimeStep(lambda, order) + default_arg('order', 4); + switch order + case 4 + k = time.rk4.get_rk4_time_step(lambda); + otherwise + error('Time-step function not available for this order'); + end + end + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+time/ExplicitRungeKuttaSecondOrderDiscreteData.m Sun Dec 23 14:39:31 2018 +0100 @@ -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/.hgtags Wed Dec 12 23:16:44 2018 +0100 +++ b/.hgtags Sun Dec 23 14:39:31 2018 +0100 @@ -2,3 +2,4 @@ 0776fa4754ff0c1918f6e1278c66f48c62d05736 grids0.1 b723495cdb2f96314d7b3f0aa79723a7dc088c7d v0.2 08f3ffe63f484d02abce8df4df61e826f568193f elastic1.0 +08f3ffe63f484d02abce8df4df61e826f568193f Heimisson2018
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscr.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,130 @@ + +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); + +% Return zeros if x0 is outside grid +if(x_0in < x(1) || x_0in > x(end) ) + + ret = zeros(size(x)); + +else + + fnorm = diag(H); + eta = abs(x-x_0in); + tot = m_order+s_order; + S = []; + M = []; + + % Get interior grid spacing + middle = floor(m/2); + h = x(middle+1) - x(middle); + + poss = find(tot*h/2 >= eta); + + % Ensure that poss is not too long + if length(poss) == (tot + 2) + poss = poss(2:end-1); + elseif length(poss) == (tot + 1) + poss = poss(1:end-1); + end + + % Use first tot grid points + if length(poss)<tot && x_0in < x(1) + ceil(tot/2)*h; + index=1:tot; + pol=(x(1:tot)-x(1))/(x(tot)-x(1)); + x_0=(x_0in-x(1))/(x(tot)-x(1)); + norm=fnorm(1:tot)/h; + + % Use last tot grid points + elseif length(poss)<tot && x_0in > x(end) - ceil(tot/2)*h; + index = length(x)-tot+1:length(x); + pol = (x(end-tot+1:end)-x(end-tot+1))/(x(end)-x(end-tot+1)); + norm = fnorm(end-tot+1:end)/h; + x_0 = (x_0in-x(end-tot+1))/(x(end)-x(end-tot+1)); + + % Interior, compensate for round-off errors. + elseif length(poss) < tot + if poss(end)<m + poss = [poss; poss(end)+1]; + else + poss = [poss(1)-1; poss]; + end + 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; + index = poss; + + % Interior + 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; + index = poss; + end + + h_pol = pol(2)-pol(1); + b = zeros(m_order+s_order,1); + + for i = 1:m_order + b(i,1) = x_0^(i-1); + end + + for i = 1:(m_order+s_order) + for j = 1:m_order + M(j,i) = pol(i)^(j-1)*h_pol*norm(i); + end + end + + for i = 1:(m_order+s_order) + for j = 1:s_order + S(j,i) = (-1)^(i-1)*pol(i)^(j-1); + end + end + + A = [M;S]; + + d = A\b; + ret = x*0; + ret(index) = d/h*h_pol; +end + +end + + + + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/diracDiscrTest.m Sun Dec 23 14:39:31 2018 +0100 @@ -0,0 +1,595 @@ +function tests = diracDiscrTest() + tests = functiontests(localfunctions); +end + +function testLeftGP(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] = setup1D(order, mom_cond); + + % Test left boundary grid points + x0s = xl + [0, h, 2*h]; + + 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 testLeftRandom(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] = setup1D(order, mom_cond); + + % Test random points near left boundary + x0s = xl + 2*h*rand(1,10); + + 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 testRightGP(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] = setup1D(order, mom_cond); + + % Test right boundary grid points + x0s = xr-[0, h, 2*h]; + + 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 testRightRandom(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] = setup1D(order, mom_cond); + + % Test random points near right boundary + x0s = xr - 2*h*rand(1,10); + + 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 testInteriorGP(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] = setup1D(order, mom_cond); + + % Test interior grid points + m_half = round(m/2); + x0s = xl + (m_half-1:m_half+1)*h; + + 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 testInteriorRandom(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] = setup1D(order, mom_cond); + + % Test random points in interior + x0s = (xl+2*h) + (xr-xl-4*h)*rand(1,20); + + 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 + +% x0 outside grid should yield 0 integral! +function testX0OutsideGrid(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] = setup1D(order, mom_cond); + + % Test points outisde grid + x0s = [xl-1.1*h, xr+1.1*h]; + + 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 - 0); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllGP(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] = setup1D(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 testHalfGP(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] = setup1D(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 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); + [xlims, ylims, m, x, X, ~, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + + % Test all grid points + x0s = X; + + for j = 1:length(fs) + f = fs{j}; + fx = f(X(:,1), X(:,2)); + for i = 1:length(x0s) + x0 = x0s(i,:); + delta = diracDiscr(x0, x, mom_cond, 0, H); + integral = delta'*H_global*fx; + err = abs(integral - f(x0(1), x0(2))); + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +function testAllRandom2D(testCase) + + orders = [2, 4, 6]; + mom_conds = orders; + + for o = 1:length(orders) + order = orders(o); + mom_cond = mom_conds(o); + [xlims, ylims, m, x, X, h, H, fs] = setup2D(order, mom_cond); + H_global = kron(H{1}, H{2}); + + 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(:,1), X(:,2)); + 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 + err = abs(integral - 0); + else + err = abs(integral - f(x0(1), x0(2))); + end + testCase.verifyLessThan(err, 1e-12); + end + end + end +end + +%=============== 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); + [xlims, ylims, zlims, m, x, X, h, H, fs] = setup3D(order, mom_cond); + H_global = kron(kron(H{1}, H{2}), H{3}); + + % Test all grid points + x0s = X; + + 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; + 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 + end +end + + +% ====================================================== +% ============== Setup functions ======================= +% ====================================================== +function [xl, xr, m, h, x, H, fs] = setup1D(order, mom_cond) + + % Grid + xl = -3; + xr = 900; + L = xr-xl; + m = 101; + h = (xr-xl)/(m-1); + g = grid.equidistant(m, {xl, xr}); + x = g.points(); + + % Quadrature + ops = sbp.D2Standard(m, {xl, xr}, order); + H = ops.H; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x) (x/L).^p; + end + +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 + xl = -3; + xr = 900; + L = xr-xl; + m = 101; + [~, g_dual] = grid.primalDual1D(m, {xl, xr}); + x = g_dual.points(); + h = g_dual.h; + + % Quadrature + ops = sbp.D1Staggered(m, {xl, xr}, order); + H = ops.H_dual; + + % Moment conditions + fs = cell(mom_cond,1); + for p = 0:mom_cond-1 + fs{p+1} = @(x) (x/L).^p; + end + +end