Mercurial > repos > public > sbplib
changeset 1227:02dfe3a56742 feature/poroelastic
Add Upwind ElasticAnisotropic schemes. Seem to work really well!
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Sat, 16 Nov 2019 14:26:06 -0800 |
| parents | e1d4cb8b5309 |
| children | 8fc2f9a4c882 |
| files | +scheme/Elastic2dCurvilinearAnisotropicUpwind.m +scheme/Elastic2dVariableAnisotropicUpwind.m |
| diffstat | 2 files changed, 1104 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Elastic2dCurvilinearAnisotropicUpwind.m Sat Nov 16 14:26:06 2019 -0800 @@ -0,0 +1,462 @@ +classdef Elastic2dCurvilinearAnisotropicUpwind < scheme.Scheme + +% Discretizes the elastic wave equation: +% rho u_{i,tt} = dj C_{ijkl} dk u_j +% in curvilinear coordinates. +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. +% Assumes fully compatible operators. + + 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 variable coefficients + J, Ji + RHO % Density + C % Elastic stiffness tensor + + D % Total operator + + Dx, Dy % Physical derivatives + sigma % Cell matrix of physical stress operators + n_w, n_e, n_s, n_n % Physical normals + + % Boundary operators in cell format, used for BC + T_w, T_e, T_s, T_n + + % Traction operators + tau_w, tau_e, tau_s, tau_n % Return vector field + tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field + tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field + + % Inner products + H + + % Boundary inner products (for scalar field) + H_w, H_e, H_s, H_n + + % Surface Jacobian vectors + s_w, s_e, s_s, s_n + + % Boundary restriction operators + e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary + e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary + e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary + e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field + en_w, en_e, en_s, en_n % Act on vector field, return normal component + + % E{i}^T picks out component i + E + + % Elastic2dVariableAnisotropic object for reference domain + refObj + end + + methods + + % The coefficients can either be function handles or grid functions + % optFlag -- if true, extra computations are performed, which may be helpful for optimization. + function obj = Elastic2dCurvilinearAnisotropicUpwind(g, order, rho, C, opSet, optFlag) + default_arg('rho', @(x,y) 0*x+1); + default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible}); + default_arg('optFlag', false); + dim = 2; + + C_default = cell(dim,dim,dim,dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + C_default{i,j,k,l} = @(x,y) 0*x + 1; + end + end + end + end + default_arg('C', C_default); + + assert(isa(g, 'grid.Curvilinear')); + + if isa(rho, 'function_handle') + rho = grid.evalOn(g, rho); + end + + C_mat = cell(dim,dim,dim,dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + if isa(C{i,j,k,l}, 'function_handle') + C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l}); + end + C_mat{i,j,k,l} = spdiag(C{i,j,k,l}); + end + end + end + end + obj.C = C_mat; + + m = g.size(); + m_tot = g.N(); + + % 1D operators + m_u = m(1); + m_v = m(2); + ops_u = opSet{1}(m_u, {0, 1}, order); + ops_v = opSet{2}(m_v, {0, 1}, order); + + h_u = ops_u.h; + h_v = ops_v.h; + + I_u = speye(m_u); + I_v = speye(m_v); + + D1_u = ops_u.D1; + H_u = ops_u.H; + Hi_u = ops_u.HI; + e_l_u = ops_u.e_l; + e_r_u = ops_u.e_r; + d1_l_u = ops_u.d1_l; + d1_r_u = ops_u.d1_r; + + D1_v = ops_v.D1; + H_v = ops_v.H; + Hi_v = ops_v.HI; + e_l_v = ops_v.e_l; + e_r_v = ops_v.e_r; + d1_l_v = ops_v.d1_l; + d1_r_v = ops_v.d1_r; + + + % Logical operators + Du = kr(D1_u,I_v); + Dv = kr(I_u,D1_v); + + e_w = kr(e_l_u,I_v); + e_e = kr(e_r_u,I_v); + e_s = kr(I_u,e_l_v); + e_n = kr(I_u,e_r_v); + + % Metric coefficients + coords = g.points(); + x = coords(:,1); + y = coords(:,2); + + x_u = Du*x; + x_v = Dv*x; + y_u = Du*y; + y_v = Dv*y; + + J = x_u.*y_v - x_v.*y_u; + + K = cell(dim, dim); + K{1,1} = y_v./J; + K{1,2} = -y_u./J; + K{2,1} = -x_v./J; + K{2,2} = x_u./J; + + % Physical derivatives + obj.Dx = spdiag( y_v./J)*Du + spdiag(-y_u./J)*Dv; + obj.Dy = spdiag(-x_v./J)*Du + spdiag( x_u./J)*Dv; + + % Wrap around Aniosotropic Cartesian + rho_tilde = J.*rho; + + PHI = cell(dim,dim,dim,dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + PHI{i,j,k,l} = 0*C{i,j,k,l}; + for m = 1:dim + for n = 1:dim + PHI{i,j,k,l} = PHI{i,j,k,l} + J.*K{m,i}.*C{m,j,n,l}.*K{n,k}; + end + end + end + end + end + end + + gRef = grid.equidistant([m_u, m_v], {0,1}, {0,1}); + refObj = scheme.Elastic2dVariableAnisotropicUpwind(gRef, order, rho_tilde, PHI, []); + + %---- Set object properties ------ + obj.RHO = spdiag(rho); + + % Volume quadrature + obj.J = spdiag(J); + obj.Ji = spdiag(1./J); + obj.H = obj.J*kr(H_u,H_v); + + % Boundary quadratures + s_w = sqrt((e_w'*x_v).^2 + (e_w'*y_v).^2); + s_e = sqrt((e_e'*x_v).^2 + (e_e'*y_v).^2); + s_s = sqrt((e_s'*x_u).^2 + (e_s'*y_u).^2); + s_n = sqrt((e_n'*x_u).^2 + (e_n'*y_u).^2); + obj.s_w = s_w; + obj.s_e = s_e; + obj.s_s = s_s; + obj.s_n = s_n; + + obj.H_w = H_v*spdiag(s_w); + obj.H_e = H_v*spdiag(s_e); + obj.H_s = H_u*spdiag(s_s); + obj.H_n = H_u*spdiag(s_n); + + % Restriction operators + obj.e_w = refObj.e_w; + obj.e_e = refObj.e_e; + obj.e_s = refObj.e_s; + obj.e_n = refObj.e_n; + + % Adapt things from reference object + obj.D = refObj.D; + obj.E = refObj.E; + + obj.e1_w = refObj.e1_w; + obj.e1_e = refObj.e1_e; + obj.e1_s = refObj.e1_s; + obj.e1_n = refObj.e1_n; + + obj.e2_w = refObj.e2_w; + obj.e2_e = refObj.e2_e; + obj.e2_s = refObj.e2_s; + obj.e2_n = refObj.e2_n; + + obj.e_scalar_w = refObj.e_scalar_w; + obj.e_scalar_e = refObj.e_scalar_e; + obj.e_scalar_s = refObj.e_scalar_s; + obj.e_scalar_n = refObj.e_scalar_n; + + e1_w = obj.e1_w; + e1_e = obj.e1_e; + e1_s = obj.e1_s; + e1_n = obj.e1_n; + + e2_w = obj.e2_w; + e2_e = obj.e2_e; + e2_s = obj.e2_s; + e2_n = obj.e2_n; + + obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')'; + obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')'; + obj.tau1_s = (spdiag(1./s_s)*refObj.tau1_s')'; + obj.tau1_n = (spdiag(1./s_n)*refObj.tau1_n')'; + + obj.tau2_w = (spdiag(1./s_w)*refObj.tau2_w')'; + obj.tau2_e = (spdiag(1./s_e)*refObj.tau2_e')'; + obj.tau2_s = (spdiag(1./s_s)*refObj.tau2_s')'; + obj.tau2_n = (spdiag(1./s_n)*refObj.tau2_n')'; + + obj.tau_w = (refObj.e_w'*obj.e1_w*obj.tau1_w')' + (refObj.e_w'*obj.e2_w*obj.tau2_w')'; + obj.tau_e = (refObj.e_e'*obj.e1_e*obj.tau1_e')' + (refObj.e_e'*obj.e2_e*obj.tau2_e')'; + obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')'; + obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')'; + + % Physical normals + e_w = obj.e_scalar_w; + e_e = obj.e_scalar_e; + e_s = obj.e_scalar_s; + e_n = obj.e_scalar_n; + + e_w_vec = obj.e_w; + e_e_vec = obj.e_e; + e_s_vec = obj.e_s; + e_n_vec = obj.e_n; + + nu_w = [-1,0]; + nu_e = [1,0]; + nu_s = [0,-1]; + nu_n = [0,1]; + + obj.n_w = cell(2,1); + obj.n_e = cell(2,1); + obj.n_s = cell(2,1); + obj.n_n = cell(2,1); + + n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2))); + n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2))); + obj.n_w{1} = spdiag(n_w_1); + obj.n_w{2} = spdiag(n_w_2); + + n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2))); + n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2))); + obj.n_e{1} = spdiag(n_e_1); + obj.n_e{2} = spdiag(n_e_2); + + n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2))); + n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2))); + obj.n_s{1} = spdiag(n_s_1); + obj.n_s{2} = spdiag(n_s_2); + + n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2))); + n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2))); + obj.n_n{1} = spdiag(n_n_1); + obj.n_n{2} = spdiag(n_n_2); + + % Operators that extract the normal component + obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')'; + obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')'; + obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')'; + obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')'; + + % Stress operators + sigma = cell(dim, dim); + D1 = {obj.Dx, obj.Dy}; + E = obj.E; + N = length(obj.RHO); + for i = 1:dim + for j = 1:dim + sigma{i,j} = sparse(N,2*N); + for k = 1:dim + for l = 1:dim + sigma{i,j} = sigma{i,j} + obj.C{i,j,k,l}*D1{k}*E{l}'; + end + end + end + end + obj.sigma = sigma; + + % Misc. + obj.refObj = refObj; + obj.m = refObj.m; + obj.h = refObj.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. Can also be e.g. + % {'normal', 'd'} or {'tangential', 't'} for conditions on + % tangential/normal 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. + + % For displacement bc: + % bc = {comp, 'd', dComps}, + % where + % dComps = vector of components with displacement BC. Default: 1:dim. + % In this way, we can specify one BC at a time even though the SATs depend on all BC. + function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) + default_arg('tuning', 1.0); + assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); + + [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning); + + type = bc{2}; + + switch type + case {'F','f','Free','free','traction','Traction','t','T'} + s = obj.(['s_' boundary]); + s = spdiag(s); + penalty = penalty*s; + end + end + + % type Struct that specifies the interface coupling. + % Fields: + % -- tuning: penalty strength, defaults to 1.0 + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.tuning = 1.0; + defaultType.interpolation = 'none'; + default_struct('type', defaultType); + + [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type); + end + + % Returns h11 for the boundary specified by the string boundary. + % op -- string + function h11 = getBorrowing(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + switch boundary + case {'w','e'} + h11 = obj.refObj.h11{1}; + case {'s', 'n'} + h11 = obj.refObj.h11{2}; + end + end + + % Returns the outward unit normal vector for the boundary specified by the string boundary. + % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2. + function n = getNormal(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + n = obj.(['n_' boundary]); + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'}) + + o = obj.([op, '_', boundary]); + + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + function o = getBoundaryOperatorForScalarField(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e'}) + + switch op + + case 'e' + o = obj.(['e_scalar', '_', boundary]); + end + + end + + % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary. + % Formula: tau_i = T_ij u_j + % op -- string + function T = getBoundaryTractionOperator(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + T = obj.(['T', '_', boundary]); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary unknowns + % + % boundary -- string + function H = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H = obj.getBoundaryQuadratureForScalarField(boundary); + I_dim = speye(obj.dim, obj.dim); + H = kron(H, I_dim); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary grid points + % + % boundary -- string + function H_b = getBoundaryQuadratureForScalarField(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H_b = obj.(['H_', boundary]); + 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/Elastic2dVariableAnisotropicUpwind.m Sat Nov 16 14:26:06 2019 -0800 @@ -0,0 +1,642 @@ +classdef Elastic2dVariableAnisotropicUpwind < scheme.Scheme + +% Discretizes the elastic wave equation: +% rho u_{i,tt} = dj C_{ijkl} dk u_j +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. +% Assumes fully compatible operators + + 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 variable coefficients + RHO, RHOi, RHOi_kron % Density + C % Elastic stiffness tensor + + D % Total operator + Dp, Dm % First derivatives + + % Boundary operators in cell format, used for BC + T_w, T_e, T_s, T_n + + % Traction operators + tau_w, tau_e, tau_s, tau_n % Return vector field + tau1_w, tau1_e, tau1_s, tau1_n % Return scalar field + tau2_w, tau2_e, tau2_s, tau2_n % Return scalar field + + % Inner products + H, Hi, Hi_kron, H_1D + + % Boundary inner products (for scalar field) + H_w, H_e, H_s, H_n + + % Boundary restriction operators + e_w, e_e, e_s, e_n % Act on vector field, return vector field at boundary + e1_w, e1_e, e1_s, e1_n % Act on vector field, return scalar field at boundary + e2_w, e2_e, e2_s, e2_n % Act on vector field, return scalar field at boundary + e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field + + % E{i}^T picks out component i + E + + % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant. + h11 % First entry in norm matrix + + end + + methods + + % The coefficients can either be function handles or grid functions + % optFlag -- if true, extra computations are performed, which may be helpful for optimization. + function obj = Elastic2dVariableAnisotropicUpwind(g, order, rho, C, opSet, optFlag) + default_arg('rho', @(x,y) 0*x+1); + default_arg('opSet',{@sbp.D1Upwind, @sbp.D1Upwind}); + default_arg('optFlag', false); + dim = 2; + + C_default = cell(dim,dim,dim,dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + C_default{i,j,k,l} = @(x,y) 0*x + 1; + end + end + end + end + default_arg('C', C_default); + assert(isa(g, 'grid.Cartesian')) + + if isa(rho, 'function_handle') + rho = grid.evalOn(g, rho); + end + + C_mat = cell(dim,dim,dim,dim); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + if isa(C{i,j,k,l}, 'function_handle') + C{i,j,k,l} = grid.evalOn(g, C{i,j,k,l}); + end + C_mat{i,j,k,l} = spdiag(C{i,j,k,l}); + end + end + end + end + obj.C = C_mat; + + m = g.size(); + m_tot = g.N(); + 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); + h = zeros(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), lim{i}, order); + h(i) = ops{i}.h; + end + + % Borrowing constants + for i = 1:dim + obj.h11{i} = ops{i}.H(1,1); + end + + I = cell(dim,1); + Dp = cell(dim,1); + Dm = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_0 = cell(dim,1); + e_m = cell(dim,1); + d1_0 = cell(dim,1); + d1_m = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + Dp{i} = ops{i}.Dp; + Dm{i} = ops{i}.Dm; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_0{i} = ops{i}.e_l; + e_m{i} = ops{i}.e_r; + d1_0{i} = (ops{i}.e_l' * Dm{i})'; + d1_m{i} = (ops{i}.e_r' * Dm{i})'; + end + + %====== Assemble full operators ======== + I_dim = speye(dim, dim); + RHO = spdiag(rho); + obj.RHO = RHO; + obj.RHOi = inv(RHO); + obj.RHOi_kron = kron(obj.RHOi, I_dim); + + obj.Dp = cell(dim,1); + obj.Dm = cell(dim,1); + + % D1 + obj.Dp{1} = kron(Dp{1},I{2}); + obj.Dp{2} = kron(I{1},Dp{2}); + obj.Dm{1} = kron(Dm{1},I{2}); + obj.Dm{2} = kron(I{1},Dm{2}); + + % Boundary restriction operators + e_l = cell(dim,1); + e_r = cell(dim,1); + e_l{1} = kron(e_0{1}, I{2}); + e_l{2} = kron(I{1}, e_0{2}); + e_r{1} = kron(e_m{1}, I{2}); + e_r{2} = kron(I{1}, e_m{2}); + + e_scalar_w = e_l{1}; + e_scalar_e = e_r{1}; + e_scalar_s = e_l{2}; + e_scalar_n = e_r{2}; + + e_w = kron(e_scalar_w, I_dim); + e_e = kron(e_scalar_e, I_dim); + e_s = kron(e_scalar_s, I_dim); + e_n = kron(e_scalar_n, I_dim); + + % 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; + + e1_w = (e_scalar_w'*E{1}')'; + e1_e = (e_scalar_e'*E{1}')'; + e1_s = (e_scalar_s'*E{1}')'; + e1_n = (e_scalar_n'*E{1}')'; + + e2_w = (e_scalar_w'*E{2}')'; + e2_e = (e_scalar_e'*E{2}')'; + e2_s = (e_scalar_s'*E{2}')'; + e2_n = (e_scalar_n'*E{2}')'; + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_w = H{2}; + obj.H_e = H{2}; + obj.H_s = H{1}; + obj.H_n = H{1}; + obj.H_1D = {H{1}, H{2}}; + + % Differentiation matrix D (without SAT) + Dp = obj.Dp; + Dm = obj.Dm; + D = sparse(dim*m_tot,dim*m_tot); + for i = 1:dim + for j = 1:dim + for k = 1:dim + for l = 1:dim + D = D + E{j}*Dp{i}*C_mat{i,j,k,l}*Dm{k}*E{l}'; + end + end + end + end + D = obj.RHOi_kron*D; + obj.D = D; + %=========================================%' + + % Numerical traction operators for BC. + % + % Formula at boundary j: % tau^{j}_i = sum_l T^{j}_{il} u_l + % + T_l = cell(dim,1); + T_r = cell(dim,1); + tau_l = cell(dim,1); + tau_r = cell(dim,1); + + D1 = obj.Dm; + + % Boundary j + 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); + + [~, n_l] = size(e_l{j}); + [~, n_r] = size(e_r{j}); + + % Traction component i + for i = 1:dim + tau_l{j}{i} = sparse(dim*m_tot, n_l); + tau_r{j}{i} = sparse(dim*m_tot, n_r); + + % Displacement component l + for l = 1:dim + T_l{j}{i,l} = sparse(m_tot, n_l); + T_r{j}{i,l} = sparse(m_tot, n_r); + + % Derivative direction k + for k = 1:dim + T_l{j}{i,l} = T_l{j}{i,l} ... + - (e_l{j}'*C_mat{j,i,k,l}*D1{k})'; + T_r{j}{i,l} = T_r{j}{i,l} ... + + (e_r{j}'*C_mat{j,i,k,l}*D1{k})'; + end + tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')'; + tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')'; + end + end + end + + % Traction tensors, T_ij + obj.T_w = T_l{1}; + obj.T_e = T_r{1}; + obj.T_s = T_l{2}; + obj.T_n = T_r{2}; + + % Restriction operators + obj.e_w = e_w; + obj.e_e = e_e; + obj.e_s = e_s; + obj.e_n = e_n; + + obj.e1_w = e1_w; + obj.e1_e = e1_e; + obj.e1_s = e1_s; + obj.e1_n = e1_n; + + obj.e2_w = e2_w; + obj.e2_e = e2_e; + obj.e2_s = e2_s; + obj.e2_n = e2_n; + + obj.e_scalar_w = e_scalar_w; + obj.e_scalar_e = e_scalar_e; + obj.e_scalar_s = e_scalar_s; + obj.e_scalar_n = e_scalar_n; + + % First component of traction + obj.tau1_w = tau_l{1}{1}; + obj.tau1_e = tau_r{1}{1}; + obj.tau1_s = tau_l{2}{1}; + obj.tau1_n = tau_r{2}{1}; + + % Second component of traction + obj.tau2_w = tau_l{1}{2}; + obj.tau2_e = tau_r{1}{2}; + obj.tau2_s = tau_l{2}{2}; + obj.tau2_n = tau_r{2}{2}; + + % Traction vectors + obj.tau_w = (e_w'*e1_w*obj.tau1_w')' + (e_w'*e2_w*obj.tau2_w')'; + obj.tau_e = (e_e'*e1_e*obj.tau1_e')' + (e_e'*e2_e*obj.tau2_e')'; + obj.tau_s = (e_s'*e1_s*obj.tau1_s')' + (e_s'*e2_s*obj.tau2_s')'; + obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')'; + + % Kroneckered norms and coefficients + 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. Can also be e.g. + % {'normal', 'd'} or {'tangential', 't'} for conditions on + % tangential/normal 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. + + % For displacement bc: + % bc = {comp, 'd', dComps}, + % where + % dComps = vector of components with displacement BC. Default: 1:dim. + % In this way, we can specify one BC at a time even though the SATs depend on all BC. + function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning) + default_arg('tuning', 1.0); + + assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' ); + comp = bc{1}; + type = bc{2}; + if ischar(comp) + comp = obj.getComponent(comp, boundary); + end + + e = obj.getBoundaryOperatorForScalarField('e', boundary); + tau = obj.getBoundaryOperator(['tau' num2str(comp)], boundary); + T = obj.getBoundaryTractionOperator(boundary); + h11 = obj.getBorrowing(boundary); + H_gamma = obj.getBoundaryQuadratureForScalarField(boundary); + nu = obj.getNormal(boundary); + + E = obj.E; + Hi = obj.Hi; + RHOi = obj.RHOi; + C = obj.C; + + dim = obj.dim; + m_tot = obj.grid.N(); + + % Preallocate + [~, col] = size(tau); + closure = sparse(dim*m_tot, dim*m_tot); + penalty = sparse(dim*m_tot, col); + + j = comp; + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet','displacement','Displacement'} + + if numel(bc) >= 3 + dComps = bc{3}; + else + dComps = 1:dim; + end + + % Loops over components that Dirichlet penalties end up on + % Y: symmetrizing part of penalty + % Z: symmetric part of penalty + % X = Y + Z. + + % Nonsymmetric part goes on all components to + % yield traction in discrete energy rate + for i = 1:dim + Y = T{j,i}'; + X = e*Y; + closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); + penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; + end + + % Symmetric part only required on components with displacement BC. + % (Otherwise it's not symmetric.) + for i = dComps + Z = sparse(m_tot, m_tot); + for l = 1:dim + for k = 1:dim + Z = Z + nu(l)*C{l,i,k,j}*nu(k); + end + end + Z = -tuning*dim/h11*Z; + X = Z; + closure = closure + E{i}*RHOi*Hi*X'*e*H_gamma*(e'*E{j}' ); + penalty = penalty - E{i}*RHOi*Hi*X'*e*H_gamma; + end + + % Free boundary condition + case {'F','f','Free','free','traction','Traction','t','T'} + closure = closure - E{j}*RHOi*Hi*e*H_gamma*tau'; + penalty = penalty + E{j}*RHOi*Hi*e*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + % type Struct that specifies the interface coupling. + % Fields: + % -- tuning: penalty strength, defaults to 1.0 + % -- interpolation: type of interpolation, default 'none' + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary,type) + + defaultType.tuning = 1.0; + 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 + + u = obj; + v = neighbour_scheme; + + % Operators, u side + e_u = u.getBoundaryOperatorForScalarField('e', boundary); + tau_u = u.getBoundaryOperator('tau', boundary); + h11_u = u.getBorrowing(boundary); + nu_u = u.getNormal(boundary); + + E_u = u.E; + C_u = u.C; + m_tot_u = u.grid.N(); + + % Operators, v side + e_v = v.getBoundaryOperatorForScalarField('e', neighbour_boundary); + tau_v = v.getBoundaryOperator('tau', neighbour_boundary); + h11_v = v.getBorrowing(neighbour_boundary); + nu_v = v.getNormal(neighbour_boundary); + + E_v = v.E; + C_v = v.C; + m_tot_v = v.grid.N(); + + % Operators that are only required for own domain + Hi = u.Hi_kron; + RHOi = u.RHOi_kron; + e_kron = u.getBoundaryOperator('e', boundary); + T_u = u.getBoundaryTractionOperator(boundary); + + % Shared operators + H_gamma = u.getBoundaryQuadratureForScalarField(boundary); + H_gamma_kron = u.getBoundaryQuadrature(boundary); + dim = u.dim; + + % Preallocate + [~, m_int] = size(H_gamma); + closure = sparse(dim*m_tot_u, dim*m_tot_u); + penalty = sparse(dim*m_tot_u, dim*m_tot_v); + + % ---- Continuity of displacement ------ + + % Y: symmetrizing part of penalty + % Z: symmetric part of penalty + % X = Y + Z. + + % Loop over components to couple across interface + for j = 1:dim + + % Loop over components that penalties end up on + for i = 1:dim + Y = 1/2*T_u{j,i}'; + Z_u = sparse(m_int, m_int); + Z_v = sparse(m_int, m_int); + for l = 1:dim + for k = 1:dim + Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u; + Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v; + end + end + Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v ); + X = Y + Z*e_u'; + closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}'; + penalty = penalty - E_u{i}*X'*H_gamma*e_v'*E_v{j}'; + end + end + + % ---- Continuity of traction ------ + closure = closure - 1/2*e_kron*H_gamma_kron*tau_u'; + penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v'; + + % ---- Multiply by inverse of density x quadraure ---- + closure = RHOi*Hi*closure; + penalty = RHOi*Hi*penalty; + + end + + function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type) + error('Non-conforming interfaces not implemented yet.'); + end + + % Returns the component number that is the tangential/normal component + % at the specified boundary + function comp = getComponent(obj, comp_str, boundary) + assertIsMember(comp_str, {'normal', 'tangential'}); + assertIsMember(boundary, {'w', 'e', 's', 'n'}); + + switch boundary + case {'w', 'e'} + switch comp_str + case 'normal' + comp = 1; + case 'tangential' + comp = 2; + end + case {'s', 'n'} + switch comp_str + case 'normal' + comp = 2; + case 'tangential' + comp = 1; + end + end + end + + % Returns h11 for the boundary specified by the string boundary. + % op -- string + function h11 = getBorrowing(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + switch boundary + case {'w','e'} + h11 = obj.h11{1}; + case {'s', 'n'} + h11 = obj.h11{2}; + end + end + + % Returns the outward unit normal vector for the boundary specified by the string boundary. + function nu = getNormal(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + switch boundary + case 'w' + nu = [-1,0]; + case 'e' + nu = [1,0]; + case 's' + nu = [0,-1]; + case 'n' + nu = [0,1]; + end + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + function o = getBoundaryOperator(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}) + + switch op + case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'} + o = obj.([op, '_', boundary]); + end + + end + + % Returns the boundary operator op for the boundary specified by the string boundary. + % op -- string + function o = getBoundaryOperatorForScalarField(obj, op, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + assertIsMember(op, {'e'}) + + switch op + + case 'e' + o = obj.(['e_scalar', '_', boundary]); + end + + end + + % Returns the boundary operator T_ij (cell format) for the boundary specified by the string boundary. + % Formula: tau_i = T_ij u_j + % op -- string + function T = getBoundaryTractionOperator(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + T = obj.(['T', '_', boundary]); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary unknowns + % + % boundary -- string + function H = getBoundaryQuadrature(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H = obj.getBoundaryQuadratureForScalarField(boundary); + I_dim = speye(obj.dim, obj.dim); + H = kron(H, I_dim); + end + + % Returns square boundary quadrature matrix, of dimension + % corresponding to the number of boundary grid points + % + % boundary -- string + function H_b = getBoundaryQuadratureForScalarField(obj, boundary) + assertIsMember(boundary, {'w', 'e', 's', 'n'}) + + H_b = obj.(['H_', boundary]); + end + + function N = size(obj) + N = obj.dim*prod(obj.m); + end + end +end