--- a/+scheme/Elastic2dVariable.m	Mon Jan 07 17:07:35 2019 +0100
+++ b/+scheme/Elastic2dVariable.m	Mon Jan 07 17:10:06 2019 +0100
@@ -30,18 +30,10 @@
         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
 
-        % 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
 
@@ -50,22 +42,38 @@
         % 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 +95,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);
@@ -183,6 +185,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);
@@ -213,7 +216,7 @@
                 end
             end
             obj.D = D;
-            %=========================================%
+            %=========================================%'
 
             % Numerical traction operators for BC.
             % Because d1 =/= e0^T*D1, the numerical tractions are different
@@ -237,20 +240,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 +269,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 +299,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 +357,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 +379,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 = -tuning*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,82 +401,59 @@
             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, tau] = obj.getBoundaryOperator({'e_tot','tau_tot'}, 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;
-            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;
+            Hi = obj.Hi_kron;
+            RHOi = obj.RHOi_kron;
 
-            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;
-            %-------------------------
+            % Penalty strength operators
+            alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha_tot', boundary);
+            alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha_tot', neighbour_boundary);
 
-            % 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;
+            closure = -RHOi*Hi*e*H_gamma*(alpha_u' + alpha_v'*e_v*e');
+            penalty = RHOi*Hi*e*H_gamma*(alpha_u'*e*e_v' + alpha_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);
+            closure = closure - 1/2*RHOi*Hi*e*H_gamma*tau';
+            penalty = penalty - 1/2*RHOi*Hi*e*H_gamma*tau_v';
 
-            % 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*RHOi*Hi*tau*H_gamma*e';
+            penalty = penalty - 1/2*RHOi*Hi*tau*H_gamma*e_v';
 
-                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};
+        end
 
-                % 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
+        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.
@@ -451,7 +478,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 +494,125 @@
                 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);
+
+                        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