sbplib: +scheme/Elastic2dVariable.m comparison

comparison +scheme/Elastic2dVariable.m @ 981:a2fcc4cf2298 feature/getBoundaryOp

Update scheme.Elastic2dVariable by copy-pasting from feature/poroelastic. This is compatible with new getBoundaryOperator in multiblock.DiffOp.

author	Martin Almquist <malmquist@stanford.edu>
date	Mon, 07 Jan 2019 17:10:06 +0100
parents	21394c78c72e
children	adbb80e60b10

comparison

equal deleted inserted replaced

-:a3accd2f1283
+:a2fcc4cf2298
 % Traction operators used for BC
 T_l, T_r
 tau_l, tau_r
-H, Hi % Inner products
+H, Hi, H_1D % 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
+% 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('lambda', @(x,y) 0*x+1);
-default_arg('mu_fun', @(x,y) 0*x+1);
+default_arg('mu', @(x,y) 0*x+1);
-default_arg('rho_fun', @(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);
+if isa(lambda, 'function_handle')
-mu = grid.evalOn(g, mu_fun);
+lambda = grid.evalOn(g, lambda);
-rho = grid.evalOn(g, rho_fun);
+end
+if isa(mu, 'function_handle')
+mu = grid.evalOn(g, mu);
+end
+if isa(rho, 'function_handle')
+rho = grid.evalOn(g, rho);
+end
 m = g.size();
 m_tot = g.N();
 h = g.scaling();
 lim = g.lim;
 ops{i} = opSet{i}(m(i), lim{i}, order);
 end
 % Borrowing constants
 for i = 1:dim
-beta = ops{i}.borrowing.R.delta_D;
+obj.theta_R{i} = h(i)*ops{i}.borrowing.R.delta_D;
-obj.H11{i} = ops{i}.borrowing.H11;
+obj.theta_H{i} = h(i)*ops{i}.borrowing.H11;
-obj.phi{i} = beta/obj.H11{i};
+obj.theta_M{i} = h(i)*ops{i}.borrowing.M.d1;
-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);
 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};
+obj.H_1D = {H{1}, H{2}};
 % E{i}^T picks out component i.
 E = cell(dim,1);
 I = speye(m_tot,m_tot);
 for i = 1:dim
 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_l{j} = cell(dim,dim);
 T_r{j} = cell(dim,dim);
 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_l{j}{i} = sparse(n_l, dim*m_tot);
-tau_r{j}{i} = sparse(m_tot,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(i,j)*LAMBDA_l*(d(i,k)*d1_l{j}' + db(i,k)*e_l{j}'*D1{k})...
--d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
+-d(j,k)*MU_l*(d(i,j)*d1_l{j}' + db(i,j)*e_l{j}'*D1{i})...
--d(i,k)*MU*e_l{j}*d1_l{j}';
+-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(i,j)*LAMBDA_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*D1{k})...
-+d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
++d(j,k)*MU_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})...
-+d(i,k)*MU*e_r{j}*d1_r{j}';
++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}';
 end
 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;
 obj.tau_r = tau_r;
 obj.m = m;
 obj.h = h;
 obj.order = order;
 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
 % Closure functions return the operators applied to the own domain to close the boundary
 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, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
 E = obj.E;
 Hi = obj.Hi;
 LAMBDA = obj.LAMBDA;
 MU = obj.MU;
 switch type
 % Dirichlet boundary condition
 case {'D','d','dirichlet','Dirichlet'}
-phi = obj.phi{j};
+alpha = obj.getBoundaryOperator('alpha', boundary);
-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};
+C = transpose(T{k,i});
-A = -d(i,k)*alpha(i,j);
+A = -tuning*e*transpose(alpha{i,k});
-B = A + C;
+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
 otherwise
 error('No such boundary condition: type = %s',type);
 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, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, boundary);
-[e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
+[e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e_tot','tau_tot'}, neighbour_boundary);
+H_gamma = obj.getBoundaryQuadrature(boundary);
 % Operators and quantities that correspond to the own domain only
-Hi = obj.Hi;
+Hi = obj.Hi_kron;
-RHOi = obj.RHOi;
+RHOi = obj.RHOi_kron;
-dim = obj.dim;
+% Penalty strength operators
-%--- Other operators ----
+alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary);
-m_tot_u = obj.grid.N();
+alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary);
-E = obj.E;
-LAMBDA_u = obj.LAMBDA;
+closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
-MU_u = obj.MU;
+penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_v');
-lambda_u = e'*LAMBDA_u*e;
-mu_u = e'*MU_u*e;
+closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
+penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
-m_tot_v = neighbour_scheme.grid.N();
-E_v = neighbour_scheme.E;
+closure = closure + 1/2*RHOi*Hi*tau*H_gamma*e';
-LAMBDA_v = neighbour_scheme.LAMBDA;
+penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
-MU_v = neighbour_scheme.MU;
-lambda_v = e_v'*LAMBDA_v*e_v;
+end
-mu_v = e_v'*MU_v*e_v;
-%-------------------------
+function [closure, penalty] = interfaceNonConforming(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+error('Non-conforming interfaces not implemented yet.');
-% 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)
 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)
+% 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'}
 j = 1;
 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
 if ~iscell(op)
 op = {op};
 end
-for i = 1:length(op)
+for k = 1:length(op)
-switch op{i}
+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);
+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) 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

Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 981:a2fcc4cf2298 feature/getBoundaryOp