sbplib: +scheme/Elastic2dVariable.m comparison

comparison +scheme/Elastic2dVariable.m @ 890:c70131daaa6e feature/d1_staggered

Merge with feature/poroelastic.

author	Martin Almquist <malmquist@stanford.edu>
date	Wed, 21 Nov 2018 18:29:29 -0800
parents	14fee299ada2
children	e30aaa4a3e09 386ef449df51

comparison

equal deleted inserted replaced

-:18e10217dca9
+:c70131daaa6e
 classdef Elastic2dVariable < scheme.Scheme
 % Discretizes the elastic wave equation:
 % rho u_{i,tt} = di lambda dj u_j + dj mu di u_j + dj mu dj u_i
 % opSet should be cell array of opSets, one per dimension. This
 % is useful if we have periodic BC in one direction.
 properties
 m % Number of points in each direction, possibly a vector
 % 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
 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;
+if isempty(lim)
+x = g.x;
+lim = cell(length(x),1);
+for i = 1:length(x)
+lim{i} = {min(x{i}), max(x{i})};
+end
+end
 % 1D operators
 ops = cell(dim,1);
 for i = 1:dim
 ops{i} = opSet{i}(m(i), lim{i}, order);
 end
 % Borrowing constants
 for i = 1:dim
-beta = ops{i}.borrowing.R.delta_D;
+obj.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;
 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);
 % Loop over components
 for i = 1:dim
 tau_l{j}{i} = sparse(m_tot,dim*m_tot);
 tau_r{j}{i} = sparse(m_tot,dim*m_tot);
 for k = 1:dim
 T_l{j}{i,k} = ...
 -d(i,j)*LAMBDA*(d(i,k)*e_l{k}*d1_l{k}' + db(i,k)*D1{k})...
 -d(j,k)*MU*(d(i,j)*e_l{i}*d1_l{i}' + db(i,j)*D1{i})...
 -d(i,k)*MU*e_l{j}*d1_l{j}';
 T_r{j}{i,k} = ...
 d(i,j)*LAMBDA*(d(i,k)*e_r{k}*d1_r{k}' + db(i,k)*D1{k})...
 +d(j,k)*MU*(d(i,j)*e_r{i}*d1_r{i}' + db(i,j)*D1{i})...
 +d(i,k)*MU*e_r{j}*d1_r{j}';
 tau_l{j}{i} = tau_l{j}{i} + T_l{j}{i,k}*E{k}';
 tau_r{j}{i} = tau_r{j}{i} + T_r{j}{i,k}*E{k}';
 end
 obj.h = h;
 obj.order = order;
 obj.grid = g;
 obj.dim = dim;
+% Used for adjoint optimization
+obj.B = cell(1,dim);
+for i = 1:dim
+obj.B{i} = zeros(m(i),m(i),m(i));
+for k = 1:m(i)
+c = sparse(m(i),1);
+c(k) = 1;
+[~, obj.B{i}(:,:,k)] = ops{i}.D2(c);
+end
+end
 end
 % Closure functions return the operators applied to the own domain to close the boundary
 % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin.
 %       boundary            is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'.
-%       type                is a cell array of strings specifying the type of boundary condition for each component.
+%       bc                  is a cell array of component and bc type, e.g. {1, 'd'} for Dirichlet condition
+%                           on the first component.
 %       data                is a function returning the data that should be applied at the boundary.
 %       neighbour_scheme    is an instance of Scheme that should be interfaced to.
 %       neighbour_boundary  is a string specifying which boundary to interface to.
-function [closure, penalty] = boundary_condition(obj, boundary, type, parameter)
+function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
-default_arg('type',{'free','free'});
+default_arg('tuning', 1.2);
-default_arg('parameter', []);
+assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
+comp = bc{1};
+type = bc{2};
 % j is the coordinate direction of the boundary
-% nj: outward unit normal component.
+j = obj.get_boundary_number(boundary);
-% nj = -1 for west, south, bottom boundaries
+[e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
-% nj = 1  for east, north, top boundaries
-[j, nj] = obj.get_boundary_number(boundary);
-switch nj
-case 1
-e = obj.e_r;
-d = obj.d1_r;
-tau = obj.tau_r{j};
-T = obj.T_r{j};
-case -1
-e = obj.e_l;
-d = obj.d1_l;
-tau = obj.tau_l{j};
-T = obj.T_l{j};
-end
 E = obj.E;
 Hi = obj.Hi;
-H_gamma = obj.H_boundary{j};
 LAMBDA = obj.LAMBDA;
 MU = obj.MU;
 RHOi = obj.RHOi;
 dim = obj.dim;
 m_tot = obj.grid.N();
-RHOi_kron = obj.RHOi_kron;
-Hi_kron = obj.Hi_kron;
 % Preallocate
 closure = sparse(dim*m_tot, dim*m_tot);
-penalty = cell(dim,1);
+penalty = sparse(dim*m_tot, m_tot/obj.m(j));
-for k = 1:dim
-penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
+k = comp;
-end
+switch type
-% Loop over components that we (potentially) have different BC on
+% Dirichlet boundary condition
-for k = 1:dim
+case {'D','d','dirichlet','Dirichlet'}
-switch type{k}
+theta_R = obj.theta_R{j};
-% Dirichlet boundary condition
+theta_H = obj.theta_H{j};
-case {'D','d','dirichlet','Dirichlet'}
+theta_M = obj.theta_M{j};
-tuning = 1.2;
+a_lambda = dim/theta_H + 1/theta_R;
-phi = obj.phi{j};
+a_mu_i = 2/theta_M;
-h = obj.h(j);
+a_mu_ij = 2/theta_H + 1/theta_R;
-h11 = obj.H11{j}*h;
-gamma = obj.gamma{j};
+d = @kroneckerDelta;  % Kronecker delta
+db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
-a_lambda = dim/h11 + 1/(h11*phi);
+alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
-a_mu_i = 2/(gamma*h);
++ d(i,j)* a_mu_i*MU ...
-a_mu_ij = 2/h11 + 1/(h11*phi);
++ db(i,j)*a_mu_ij*MU );
-d = @kroneckerDelta;  % Kronecker delta
+% Loop over components that Dirichlet penalties end up on
-db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+for i = 1:dim
-alpha = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
+C = T{k,i};
-+ d(i,j)* a_mu_i*MU ...
+A = -d(i,k)*alpha(i,j);
-+ db(i,j)*a_mu_ij*MU );
+B = A + C;
+closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
-% Loop over components that Dirichlet penalties end up on
+penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
-for i = 1:dim
+end
-C = T{k,i};
-A = -d(i,k)*alpha(i,j);
+% Free boundary condition
-B = A + C;
+case {'F','f','Free','free','traction','Traction','t','T'}
-closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' );
+closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} );
-penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma;
+penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
-end
+% Unknown boundary condition
-% Free boundary condition
+otherwise
-case {'F','f','Free','free','traction','Traction','t','T'}
+error('No such boundary condition: type = %s',type);
-closure = closure - E{k}*RHOi*Hi*e{j}*H_gamma* (e{j}'*tau{k} );
-penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma;
-% Unknown boundary condition
-otherwise
-error('No such boundary condition: type = %s',type);
-end
 end
 end
 function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
 % u denotes the solution in the own domain
 % v denotes the solution in the neighbour domain
+% Operators without subscripts are from the own domain.
 tuning = 1.2;
-% tuning = 20.2;
-error('Interface not implemented');
+% j is the coordinate direction of the boundary
+j = obj.get_boundary_number(boundary);
+j_v = neighbour_scheme.get_boundary_number(neighbour_boundary);
+% Get boundary operators
+[e, T, tau, H_gamma] = obj.get_boundary_operator({'e','T','tau','H'}, boundary);
+[e_v, tau_v] = neighbour_scheme.get_boundary_operator({'e','tau'}, neighbour_boundary);
+% Operators and quantities that correspond to the own domain only
+Hi = obj.Hi;
+RHOi = obj.RHOi;
+dim = obj.dim;
+%--- Other operators ----
+m_tot_u = obj.grid.N();
+E = obj.E;
+LAMBDA_u = obj.LAMBDA;
+MU_u = obj.MU;
+lambda_u = e'*LAMBDA_u*e;
+mu_u = e'*MU_u*e;
+m_tot_v = neighbour_scheme.grid.N();
+E_v = neighbour_scheme.E;
+LAMBDA_v = neighbour_scheme.LAMBDA;
+MU_v = neighbour_scheme.MU;
+lambda_v = e_v'*LAMBDA_v*e_v;
+mu_v = e_v'*MU_v*e_v;
+%-------------------------
+% Borrowing constants
+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};
+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_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 = tuning*(alpha_ii_u + alpha_ii_v);
+sigma_ij = tuning*(alpha_ij_u + alpha_ij_v);
+d = @kroneckerDelta;  % Kronecker delta
+db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+sigma = @(i,j) tuning*(d(i,j)*sigma_ii + db(i,j)*sigma_ij);
+% Preallocate
+closure = sparse(dim*m_tot_u, dim*m_tot_u);
+penalty = sparse(dim*m_tot_u, dim*m_tot_v);
+% Loop over components that penalties end up on
+for i = 1:dim
+closure = closure - E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e'*E{i}';
+penalty = penalty + E{i}*RHOi*Hi*e*sigma(i,j)*H_gamma*e_v'*E_v{i}';
+closure = closure - 1/2*E{i}*RHOi*Hi*e*H_gamma*e'*tau{i};
+penalty = penalty - 1/2*E{i}*RHOi*Hi*e*H_gamma*e_v'*tau_v{i};
+% Loop over components that we have interface conditions on
+for k = 1:dim
+closure = closure + 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e'*E{k}';
+penalty = penalty - 1/2*E{i}*RHOi*Hi*T{k,i}'*e*H_gamma*e_v'*E_v{k}';
+end
+end
 end
 % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
 function [j, nj] = get_boundary_number(obj, boundary)
 case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
 nj = 1;
 end
 end
-% Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+% Returns the boundary operator op for the boundary specified by the string boundary.
-function [return_op] = get_boundary_operator(obj, op, boundary)
+% op: may be a cell array of strings
+function [varargout] = get_boundary_operator(obj, op, boundary)
 switch boundary
 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
 j = 1;
 case {'s','S','south','South', 'n', 'N', 'north', 'North'}
 j = 2;
 otherwise
 error('No such boundary: boundary = %s',boundary);
 end
-switch op
+if ~iscell(op)
-case 'e'
+op = {op};
-switch boundary
+end
-case {'w','W','west','West','s','S','south','South'}
-return_op = obj.e_l{j};
+for i = 1:length(op)
-case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+switch op{i}
-return_op = obj.e_r{j};
+case 'e'
-end
+switch boundary
-case 'd'
+case {'w','W','west','West','s','S','south','South'}
-switch boundary
+varargout{i} = obj.e_l{j};
-case {'w','W','west','West','s','S','south','South'}
+case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
-return_op = obj.d1_l{j};
+varargout{i} = obj.e_r{j};
-case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+end
-return_op = obj.d1_r{j};
+case 'd'
-end
+switch boundary
-otherwise
+case {'w','W','west','West','s','S','south','South'}
-error(['No such operator: operatr = ' op]);
+varargout{i} = obj.d1_l{j};
+case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+varargout{i} = obj.d1_r{j};
+end
+case 'H'
+varargout{i} = obj.H_boundary{j};
+case 'T'
+switch boundary
+case {'w','W','west','West','s','S','south','South'}
+varargout{i} = obj.T_l{j};
+case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+varargout{i} = obj.T_r{j};
+end
+case 'tau'
+switch boundary
+case {'w','W','west','West','s','S','south','South'}
+varargout{i} = obj.tau_l{j};
+case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+varargout{i} = obj.tau_r{j};
+end
+case 'alpha'
+% alpha = alpha(i,j) is the penalty strength for displacement BC.
+tuning = 1.2;
+LAMBDA = obj.LAMBDA;
+MU = obj.MU;
+phi = obj.phi{j};
+h = obj.h(j);
+h11 = obj.H11{j}*h;
+gamma = obj.gamma{j};
+dim = obj.dim;
+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,k) d(i,k)*tuning*( d(i,j)* a_lambda*LAMBDA ...
++ d(i,j)* a_mu_i*MU ...
++ db(i,j)*a_mu_ij*MU );
+varargout{i} = alpha;
+otherwise
+error(['No such operator: operator = ' op{i}]);
+end
 end
 end
 function N = size(obj)
-N = prod(obj.m);
+N = obj.dim*prod(obj.m);
 end
 end
 end

Mercurial > repos > public > sbplib

comparison +scheme/Elastic2dVariable.m @ 890:c70131daaa6e feature/d1_staggered