--- a/+scheme/Elastic2dVariable.m	Sun Nov 04 12:36:30 2018 -0800
+++ b/+scheme/Elastic2dVariable.m	Wed Nov 21 18:29:29 2018 -0800
@@ -1,7 +1,7 @@
 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 
+% 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.
 
@@ -30,40 +30,60 @@
         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
         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('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();
 
             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);
@@ -73,10 +93,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;
+                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);
@@ -164,6 +183,7 @@
             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);
@@ -194,7 +214,7 @@
                 end
             end
             obj.D = D;
-            %=========================================%
+            %=========================================%'
 
             % Numerical traction operators for BC.
             % Because d1 =/= e0^T*D1, the numerical tractions are different
@@ -223,14 +243,14 @@
                     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} = ... 
+                        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(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} = ... 
+                        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(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}';
@@ -256,41 +276,41 @@
             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)
-            default_arg('type',{'free','free'});
-            default_arg('parameter', []);
+        function [closure, penalty] = boundary_condition(obj, boundary, bc, tuning)
+            default_arg('tuning', 1.2);
+
+            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. 
-            % nj = -1 for west, south, bottom boundaries
-            % 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
+            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;
-            H_gamma = obj.H_boundary{j};
             LAMBDA = obj.LAMBDA;
             MU = obj.MU;
             RHOi = obj.RHOi;
@@ -298,66 +318,126 @@
             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);
-            for k = 1:dim
-                penalty{k} = sparse(dim*m_tot, m_tot/obj.m(j));
-            end
+            penalty = sparse(dim*m_tot, m_tot/obj.m(j));
 
-            % Loop over components that we (potentially) have different BC on
-            for k = 1:dim
-                switch type{k}
+            k = comp;
+            switch type
+
+            % Dirichlet boundary condition
+            case {'D','d','dirichlet','Dirichlet'}
 
-                % Dirichlet boundary condition
-                case {'D','d','dirichlet','Dirichlet'}
+                theta_R = obj.theta_R{j};
+                theta_H = obj.theta_H{j};
+                theta_M = obj.theta_M{j};
 
-                    tuning = 1.2;
-                    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);
+                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 = @(i,j) tuning*( d(i,j)* a_lambda*LAMBDA ...
-                                          + d(i,j)* a_mu_i*MU ...
-                                          + db(i,j)*a_mu_ij*MU ); 
+                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};
-                        A = -d(i,k)*alpha(i,j);
-                        B = A + C;
-                        closure = closure + E{i}*RHOi*Hi*B'*e{j}*H_gamma*(e{j}'*E{k}' ); 
-                        penalty{k} = penalty{k} - E{i}*RHOi*Hi*B'*e{j}*H_gamma;
-                    end 
+                % 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;
+                    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{j}*H_gamma* (e{j}'*tau{k} ); 
-                        penalty{k} = penalty{k} + E{k}*RHOi*Hi*e{j}*H_gamma;
+            % Free boundary condition
+            case {'F','f','Free','free','traction','Traction','t','T'}
+                    closure = closure - E{k}*RHOi*Hi*e*H_gamma* (e'*tau{k} );
+                    penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
 
-                % Unknown boundary condition
-                otherwise
-                    error('No such boundary condition: type = %s',type);
-                end
+            % Unknown boundary condition
+            otherwise
+                error('No such boundary condition: type = %s',type);
             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.
@@ -380,8 +460,9 @@
             end
         end
 
-        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
-        function [return_op] = get_boundary_operator(obj, op, boundary)
+        % 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)
 
             switch boundary
                 case {'w','W','west','West', 'e', 'E', 'east', 'East'}
@@ -392,29 +473,73 @@
                     error('No such boundary: boundary = %s',boundary);
             end
 
-            switch op
-                case 'e'
-                    switch boundary
-                        case {'w','W','west','West','s','S','south','South'}
-                            return_op = obj.e_l{j};
-                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
-                            return_op = obj.e_r{j};
-                    end
-                case 'd'
-                    switch boundary
-                        case {'w','W','west','West','s','S','south','South'}
-                            return_op = obj.d1_l{j};
-                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
-                            return_op = obj.d1_r{j};
-                    end
-                otherwise
-                    error(['No such operator: operatr = ' op]);
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:length(op)
+                switch op{i}
+                    case 'e'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                varargout{i} = obj.e_l{j};
+                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                                varargout{i} = obj.e_r{j};
+                        end
+                    case 'd'
+                        switch boundary
+                            case {'w','W','west','West','s','S','south','South'}
+                                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