Mercurial > repos > public > sbplib
changeset 971:e54c2f54dbfe feature/getBoundaryOperator
Merge with feature/poroelastic. Use only the changes made to multiblock.DiffOp and scheme.Elastic2dVariable. DiffOp.getBoundaryOperator/Quadrature now use scheme methods instead of propeties.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Tue, 25 Dec 2018 07:50:07 +0100 |
| parents | 2412f407749a (current diff) 23d9ca6755be (diff) |
| children | 47b48a97c675 |
| files | +multiblock/DiffOp.m +multiblock/Grid.m +sbp/+implementations/d2_variable_2.m +sbp/+implementations/d2_variable_4.m +sbp/+implementations/d4_variable_6.m +time/+rkparameters/rk4.m +time/ExplicitRungeKuttaDiscreteData.m +time/ExplicitRungeKuttaSecondOrderDiscreteData.m .hgtags diracDiscr.m diracDiscrTest.m |
| diffstat | 2 files changed, 268 insertions(+), 102 deletions(-) [+] |
line wrap: on
line diff
--- a/+multiblock/DiffOp.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+multiblock/DiffOp.m Tue Dec 25 07:50:07 2018 +0100 @@ -129,19 +129,20 @@ % Get a boundary operator specified by opName for the given boundary/BoundaryGroup function op = getBoundaryOperator(obj, opName, boundary) + blockmatrixDiv = obj.blockmatrixDiv{1}; + switch class(boundary) case 'cell' - localOpName = [opName '_' boundary{2}]; blockId = boundary{1}; - localOp = obj.diffOps{blockId}.(localOpName); + localOp = obj.diffOps{blockId}.getBoundaryOperator(opName, boundary{2}); - div = {obj.blockmatrixDiv{1}, size(localOp,2)}; + div = {blockmatrixDiv, size(localOp,2)}; blockOp = blockmatrix.zero(div); blockOp{blockId,1} = localOp; op = blockmatrix.toMatrix(blockOp); return case 'multiblock.BoundaryGroup' - op = sparse(size(obj.D,1),0); + op = sparse(sum(blockmatrixDiv),0); for i = 1:length(boundary) op = [op, obj.getBoundaryOperator(opName, boundary{i})]; end @@ -150,14 +151,12 @@ end end + % Get square matrix that integrates the solution restricted to a boundary function op = getBoundaryQuadrature(obj, boundary) - opName = 'H'; switch class(boundary) case 'cell' - localOpName = [opName '_' boundary{2}]; blockId = boundary{1}; - op = obj.diffOps{blockId}.(localOpName); - + op = obj.diffOps{blockId}.getBoundaryQuadrature(boundary{2}); return case 'multiblock.BoundaryGroup' N = length(boundary);
--- a/+scheme/Elastic2dVariable.m Wed Dec 12 23:16:44 2018 +0100 +++ b/+scheme/Elastic2dVariable.m Tue Dec 25 07:50:07 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 @@ -295,7 +367,8 @@ % j is the coordinate direction of the boundary j = obj.get_boundary_number(boundary); - [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); + [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); + E = obj.E; Hi = obj.Hi; @@ -316,33 +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.getBoundaryOperator('alpha', boundary); % Loop over components that Dirichlet penalties end up on for i = 1:dim - C = T{k,i}; - A = -d(i,k)*alpha(i,j); - B = A + C; + C = transpose(T{k,i}); + A = -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,19 +411,40 @@ end end - function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % type Struct that specifies the interface coupling. + % Fields: + % -- tuning: penalty strength, defaults to 1.2 + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.tuning = 1.2; + defaultType.interpolation = 'none'; + default_struct('type', defaultType); + + switch type.interpolation + case {'none', ''} + [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type); + case {'op','OP'} + [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type); + otherwise + error('Unknown type of interpolation: %s ', type.interpolation); + end + end + + function [closure, penalty] = interfaceStandard(obj,boundary,neighbour_scheme,neighbour_boundary,type) + tuning = type.tuning; + % u denotes the solution in the own domain % v denotes the solution in the neighbour domain % Operators without subscripts are from the own domain. - tuning = 1.2; % j is the coordinate direction of the boundary j = obj.get_boundary_number(boundary); j_v = neighbour_scheme.get_boundary_number(neighbour_boundary); % Get boundary operators - [e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary); - [e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary); + [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary); + [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_boundary); % Operators and quantities that correspond to the own domain only Hi = obj.Hi; @@ -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) @@ -451,7 +536,8 @@ % Returns the boundary operator op for the boundary specified by the string boundary. % op: may be a cell array of strings - function [varargout] = get_boundary_operator(obj, op, boundary) + % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator() + function [varargout] = getBoundaryOperator(obj, op, boundary) switch boundary case {'w','W','west','West', 'e', 'E', 'east', 'East'} @@ -466,45 +552,126 @@ op = {op}; end - for i = 1:length(op) - switch op{i} + for k = 1:length(op) + switch op{k} case 'e' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.e_l{j}; + varargout{k} = obj.e_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.e_r{j}; + varargout{k} = obj.e_r{j}; end + + case 'e_tot' + e = obj.getBoundaryOperator('e', boundary); + I_dim = speye(obj.dim, obj.dim); + varargout{k} = kron(e, I_dim); + case 'd' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.d1_l{j}; + varargout{k} = obj.d1_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.d1_r{j}; + varargout{k} = obj.d1_r{j}; end - case 'H' - varargout{i} = obj.H_boundary{j}; + case 'T' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.T_l{j}; + varargout{k} = obj.T_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.T_r{j}; + varargout{k} = obj.T_r{j}; end + case 'tau' switch boundary case {'w','W','west','West','s','S','south','South'} - varargout{i} = obj.tau_l{j}; + varargout{k} = obj.tau_l{j}; case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} - varargout{i} = obj.tau_r{j}; + varargout{k} = obj.tau_r{j}; + end + + case 'tau_tot' + [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary); + + I_dim = speye(obj.dim, obj.dim); + e_tot = kron(e, I_dim); + E = obj.E; + tau_tot = (e_tot'*E{1}*e*tau{1}')'; + for i = 2:obj.dim + tau_tot = tau_tot + (e_tot'*E{i}*e*tau{i}')'; end + varargout{k} = tau_tot; + + case 'H' + varargout{k} = obj.H_boundary{j}; + + case 'alpha' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + e = obj.getBoundaryOperator('e', boundary); + + 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; + + case 'alpha_tot' + % alpha = alpha(i,j) is the penalty strength for displacement BC. + [e, e_tot, alpha] = obj.getBoundaryOperator({'e', 'e_tot', 'alpha'}, boundary); + E = obj.E; + [m, n] = size(alpha{1,1}); + alpha_tot = sparse(m*obj.dim, n*obj.dim); + for i = 1:obj.dim + for l = 1:obj.dim + alpha_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')'; + end + end + varargout{k} = alpha_tot; + otherwise - error(['No such operator: operator = ' op{i}]); + error(['No such operator: operator = ' op{k}]); end end end + function H = getBoundaryQuadrature(obj, boundary) + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + H = obj.H_boundary{j}; + I_dim = speye(obj.dim, obj.dim); + H = kron(H, I_dim); + end + function N = size(obj) N = obj.dim*prod(obj.m); end