--- a/+scheme/Elastic2dVariable.m	Thu Feb 17 18:55:11 2022 +0100
+++ b/+scheme/Elastic2dVariable.m	Thu Mar 10 16:54:26 2022 +0100
@@ -14,10 +14,10 @@
 
         order % Order of accuracy for the approximation
 
-        % Diagonal matrices for varible coefficients
-        LAMBDA % Variable coefficient, related to dilation
-        MU     % Shear modulus, variable coefficient
-        RHO, RHOi % Density, variable
+        % Diagonal matrices for variable coefficients
+        LAMBDA % Lame's first parameter, related to dilation
+        MU     % Shear modulus
+        RHO, RHOi, RHOi_kron % Density
 
         D % Total operator
         D1 % First derivatives
@@ -26,22 +26,28 @@
         D2_lambda
         D2_mu
 
-        % Traction operators used for BC
-        T_l, T_r
-        tau_l, tau_r
+        % Boundary operators in cell format, used for BC
+        T_w, T_e, T_s, T_n
 
-        H, Hi, H_1D % Inner products
-        e_l, e_r
+        % 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
 
-        d1_l, d1_r % Normal derivatives at the boundary
-        E % E{i}^T picks out component i
+        % Boundary inner products (for scalar field)
+        H_w, H_e, H_s, H_n
 
-        H_boundary % Boundary inner products
+        % 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
 
-        % Kroneckered norms and coefficients
-        RHOi_kron
-        Hi_kron
+        % E{i}^T picks out component i
+        E
 
         % Borrowing constants of the form gamma*h, where gamma is a dimensionless constant.
         theta_R % Borrowing (d1- D1)^2 from R
@@ -55,11 +61,13 @@
     methods
 
         % The coefficients can either be function handles or grid functions
-        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
+        % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
+        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet, optFlag)
             default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
             default_arg('lambda', @(x,y) 0*x+1);
             default_arg('mu', @(x,y) 0*x+1);
             default_arg('rho', @(x,y) 0*x+1);
+            default_arg('optFlag', false);
             dim = 2;
 
             assert(isa(g, 'grid.Cartesian'))
@@ -105,10 +113,10 @@
             D2 = cell(dim,1);
             H = cell(dim,1);
             Hi = cell(dim,1);
-            e_l = cell(dim,1);
-            e_r = cell(dim,1);
-            d1_l = cell(dim,1);
-            d1_r = 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));
@@ -116,10 +124,10 @@
                 D2{i} = ops{i}.D2;
                 H{i} =  ops{i}.H;
                 Hi{i} = ops{i}.HI;
-                e_l{i} = ops{i}.e_l;
-                e_r{i} = ops{i}.e_r;
-                d1_l{i} = ops{i}.d1_l;
-                d1_r{i} = ops{i}.d1_r;
+                e_0{i} = ops{i}.e_l;
+                e_m{i} = ops{i}.e_r;
+                d1_0{i} = ops{i}.d1_l;
+                d1_m{i} = ops{i}.d1_r;
             end
 
             %====== Assemble full operators ========
@@ -134,30 +142,64 @@
             obj.D1 = cell(dim,1);
             obj.D2_lambda = cell(dim,1);
             obj.D2_mu = cell(dim,1);
-            obj.e_l = cell(dim,1);
-            obj.e_r = cell(dim,1);
-            obj.d1_l = cell(dim,1);
-            obj.d1_r = cell(dim,1);
 
             % D1
             obj.D1{1} = kron(D1{1},I{2});
             obj.D1{2} = kron(I{1},D1{2});
 
-            % Boundary operators
-            obj.e_l{1} = kron(e_l{1},I{2});
-            obj.e_l{2} = kron(I{1},e_l{2});
-            obj.e_r{1} = kron(e_r{1},I{2});
-            obj.e_r{2} = kron(I{1},e_r{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};
+
+            I_dim = speye(dim, dim);
+            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);
 
-            obj.d1_l{1} = kron(d1_l{1},I{2});
-            obj.d1_l{2} = kron(I{1},d1_l{2});
-            obj.d1_r{1} = kron(d1_r{1},I{2});
-            obj.d1_r{2} = kron(I{1},d1_r{2});
+            % Boundary derivatives
+            d1_l = cell(dim,1);
+            d1_r = cell(dim,1);
+            d1_l{1} = kron(d1_0{1}, I{2});
+            d1_l{2} = kron(I{1}, d1_0{2});
+            d1_r{1} = kron(d1_m{1}, I{2});
+            d1_r{2} = kron(I{1}, d1_m{2});
+
+
+            % 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}')';
+
 
             % D2
             for i = 1:dim
-                obj.D2_lambda{i} = sparse(m_tot);
-                obj.D2_mu{i} = sparse(m_tot);
+                obj.D2_lambda{i} = sparse(m_tot, m_tot);
+                obj.D2_mu{i} = sparse(m_tot, m_tot);
             end
             ind = grid.funcToMatrix(g, 1:m_tot);
 
@@ -182,21 +224,12 @@
             % Quadratures
             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_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}};
 
-            % 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;
-
             % Differentiation matrix D (without SAT)
             D2_lambda = obj.D2_lambda;
             D2_mu = obj.D2_mu;
@@ -221,16 +254,14 @@
             % Numerical traction operators for BC.
             % Because d1 =/= e0^T*D1, the numerical tractions are different
             % at every boundary.
+            %
+            % Formula at boundary j: % tau^{j}_i = sum_k T^{j}_{ik} u_k
+            %
             T_l = cell(dim,1);
             T_r = cell(dim,1);
             tau_l = cell(dim,1);
             tau_r = cell(dim,1);
-            % tau^{j}_i = sum_k T^{j}_{ik} u_k
 
-            d1_l = obj.d1_l;
-            d1_r = obj.d1_r;
-            e_l = obj.e_l;
-            e_r = obj.e_r;
             D1 = obj.D1;
 
             % Loop over boundaries
@@ -250,45 +281,72 @@
 
                 % Loop over components
                 for i = 1:dim
-                    tau_l{j}{i} = sparse(n_l, dim*m_tot);
-                    tau_r{j}{i} = sparse(n_r, dim*m_tot);
+                    tau_l{j}{i} = sparse(dim*m_tot, n_l);
+                    tau_r{j}{i} = sparse(dim*m_tot, n_r);
                     for k = 1:dim
                         T_l{j}{i,k} = ...
-                        -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}';
+                        (-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_r*(d(i,k)*d1_r{j}' + db(i,k)*e_r{j}'*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_r*(d(i,j)*d1_r{j}' + db(i,j)*e_r{j}'*D1{i})...
-                        +d(i,k)*MU_r*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}';
+                        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
+            % 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.T_l = T_l;
-            obj.T_r = T_r;
-            obj.tau_l = tau_l;
-            obj.tau_r = tau_r;
+            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
-            I_dim = speye(dim);
             obj.RHOi_kron = kron(obj.RHOi, I_dim);
             obj.Hi_kron = kron(obj.Hi, I_dim);
 
@@ -300,42 +358,46 @@
             obj.dim = dim;
 
             % B, used for adjoint optimization
-            B = cell(dim, 1);
-            for i = 1:dim
-                B{i} = cell(m_tot, 1);
-            end
+            B = [];
+            if optFlag
+                B = cell(dim, 1);
+                for i = 1:dim
+                    B{i} = cell(m_tot, 1);
+                end
+
+                B0 = sparse(m_tot, m_tot);
+                for i = 1:dim
+                    for j = 1:m_tot
+                        B{i}{j} = B0;
+                    end
+                end
+
+                ind = grid.funcToMatrix(g, 1:m_tot);
 
-            for i = 1:dim
-                for j = 1:m_tot
-                    B{i}{j} = sparse(m_tot, 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
             end
-
-            ind = grid.funcToMatrix(g, 1:m_tot);
+            obj.B = B;
 
-            % 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
 
@@ -344,7 +406,9 @@
         % 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.
+        %                           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.
@@ -354,11 +418,15 @@
             assert( iscell(bc), 'The BC type must be a 2x1 cell array' );
             comp = bc{1};
             type = bc{2};
+            if ischar(comp)
+                comp = obj.getComponent(comp, boundary);
+            end
 
-            % j is the coordinate direction of the boundary
-            j = obj.get_boundary_number(boundary);
-            [e, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
-
+            e       = obj.getBoundaryOperatorForScalarField('e', boundary);
+            tau     = obj.getBoundaryOperator(['tau' num2str(comp)], boundary);
+            T       = obj.getBoundaryTractionOperator(boundary);
+            alpha   = obj.getBoundaryOperatorForScalarField('alpha', boundary);
+            H_gamma = obj.getBoundaryQuadratureForScalarField(boundary);
 
             E = obj.E;
             Hi = obj.Hi;
@@ -370,8 +438,9 @@
             m_tot = obj.grid.N();
 
             % Preallocate
+            [~, col] = size(tau);
             closure = sparse(dim*m_tot, dim*m_tot);
-            penalty = sparse(dim*m_tot, m_tot/obj.m(j));
+            penalty = sparse(dim*m_tot, col);
 
             k = comp;
             switch type
@@ -379,8 +448,6 @@
             % Dirichlet boundary condition
             case {'D','d','dirichlet','Dirichlet'}
 
-                alpha = obj.getBoundaryOperator('alpha', boundary);
-
                 % Loop over components that Dirichlet penalties end up on
                 for i = 1:dim
                     C = transpose(T{k,i});
@@ -392,7 +459,7 @@
 
             % Free boundary condition
             case {'F','f','Free','free','traction','Traction','t','T'}
-                    closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau{k}';
+                    closure = closure - E{k}*RHOi*Hi*e*H_gamma*tau';
                     penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
 
             % Unknown boundary condition
@@ -429,11 +496,11 @@
             % Operators without subscripts are from the own domain.
 
             % Get boundary operators
-            e = obj.getBoundaryOperator('e_tot', boundary);
-            tau = obj.getBoundaryOperator('tau_tot', boundary);
+            e   = obj.getBoundaryOperator('e', boundary);
+            tau = obj.getBoundaryOperator('tau', boundary);
 
-            e_v = neighbour_scheme.getBoundaryOperator('e_tot', neighbour_boundary);
-            tau_v = neighbour_scheme.getBoundaryOperator('tau_tot', neighbour_boundary);
+            e_v   = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
+            tau_v = neighbour_scheme.getBoundaryOperator('tau', neighbour_boundary);
 
             H_gamma = obj.getBoundaryQuadrature(boundary);
 
@@ -442,8 +509,8 @@
             RHOi = obj.RHOi_kron;
 
             % 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);
+            alpha_u = 1/4*tuning*obj.getBoundaryOperator('alpha', boundary);
+            alpha_v = 1/4*tuning*neighbour_scheme.getBoundaryOperator('alpha', neighbour_boundary);
 
             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');
@@ -460,103 +527,108 @@
             error('Non-conforming interfaces not implemented yet.');
         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)
-            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+        % 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'}
-                    j = 1;
-                case {'s', 'n'}
-                    j = 2;
-            end
-
-            switch boundary
-                case {'w', 's'}
-                    nj = -1;
-                case {'e', 'n'}
-                    nj = 1;
+            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 the boundary operator op for the boundary specified by the string boundary.
         % op -- string
-        % Only operators with name *_tot can be used with multiblock.DiffOp.getBoundaryOperator()
-        function [varargout] = getBoundaryOperator(obj, op, boundary)
+        function o = getBoundaryOperator(obj, op, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
-            assertIsMember(op, {'e', 'e_tot', 'd', 'T', 'tau', 'tau_tot', 'H', 'alpha', 'alpha_tot'})
-
-            switch boundary
-                case {'w', 'e'}
-                    j = 1;
-                case {'s', 'n'}
-                    j = 2;
-            end
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha', 'alpha1', 'alpha2'})
 
             switch op
-                case 'e'
-                    switch boundary
-                        case {'w', 's'}
-                            o = obj.e_l{j};
-                        case {'e', 'n'}
-                            o = obj.e_r{j};
-                    end
+
+                case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
+                    o = obj.([op, '_', boundary]);
 
-                case 'e_tot'
-                    e = obj.getBoundaryOperator('e', boundary);
-                    I_dim = speye(obj.dim, obj.dim);
-                    o = kron(e, I_dim);
+                % Yields vector-valued penalty strength given displacement BC on all components
+                case 'alpha'
+                    e               = obj.getBoundaryOperator('e', boundary);
+                    e_scalar        = obj.getBoundaryOperatorForScalarField('e', boundary);
+                    alpha_scalar    = obj.getBoundaryOperatorForScalarField('alpha', boundary);
+                    E = obj.E;
+                    [m, n] = size(alpha_scalar{1,1});
+                    alpha = sparse(m*obj.dim, n*obj.dim);
+                    for i = 1:obj.dim
+                        for l = 1:obj.dim
+                            alpha = alpha + (e'*E{i}*e_scalar*alpha_scalar{i,l}'*E{l}')';
+                        end
+                    end
+                    o = alpha;
 
-                case 'd'
-                    switch boundary
-                        case {'w', 's'}
-                            o = obj.d1_l{j};
-                        case {'e', 'n'}
-                            o = obj.d1_r{j};
-                    end
+                % Yields penalty strength for component 1 given displacement BC on all components
+                case 'alpha1'
+                    alpha   = obj.getBoundaryOperator('alpha', boundary);
+                    e       = obj.getBoundaryOperator('e', boundary);
+                    e1      = obj.getBoundaryOperator('e1', boundary);
 
-                case 'T'
-                    switch boundary
-                        case {'w', 's'}
-                            o = obj.T_l{j};
-                        case {'e', 'n'}
-                            o = obj.T_r{j};
-                    end
+                    alpha1 = (e1'*e*alpha')';
+                    o = alpha1;
+
+                % Yields penalty strength for component 2 given displacement BC on all components
+                case 'alpha2'
+                    alpha   = obj.getBoundaryOperator('alpha', boundary);
+                    e       = obj.getBoundaryOperator('e', boundary);
+                    e2      = obj.getBoundaryOperator('e2', boundary);
+
+                    alpha2 = (e2'*e*alpha')';
+                    o = alpha2;
+            end
 
-                case 'tau'
-                    switch boundary
-                        case {'w', 's'}
-                            o = obj.tau_l{j};
-                        case {'e', 'n'}
-                            o = obj.tau_r{j};
-                    end
-
-                case 'tau_tot'
-                    [e, tau] = obj.getBoundaryOperator({'e', 'tau'}, boundary);
+        end
 
-                    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
-                    o = tau_tot;
+        % 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', 'alpha'})
 
-                case 'H'
-                    o = obj.H_boundary{j};
+            switch op
+
+                case 'e'
+                    o = obj.(['e_scalar', '_', boundary]);
 
                 case 'alpha'
-                    % alpha = alpha(i,j) is the penalty strength for displacement BC.
-                    e = obj.getBoundaryOperator('e', boundary);
+
+                    % alpha{i,j} is the penalty strength on component i due to
+                    % displacement BC for component j.
+                    e = obj.getBoundaryOperatorForScalarField('e', boundary);
 
                     LAMBDA = obj.LAMBDA;
                     MU = obj.MU;
+                    dim = obj.dim;
 
-                    dim = obj.dim;
-                    theta_R = obj.theta_R{j};
-                    theta_H = obj.theta_H{j};
-                    theta_M = obj.theta_M{j};
+                    switch boundary
+                        case {'w', 'e'}
+                            k = 1;
+                        case {'s', 'n'}
+                            k = 2;
+                    end
+
+                    theta_R = obj.theta_R{k};
+                    theta_H = obj.theta_H{k};
+                    theta_M = obj.theta_M{k};
 
                     a_lambda = dim/theta_H + 1/theta_R;
                     a_mu_i = 2/theta_M;
@@ -564,53 +636,53 @@
 
                     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;
+
+                    alpha = cell(obj.dim, obj.dim);
                     for i = 1:obj.dim
-                        for l = 1:obj.dim
-                            alpha{i,l} = d(i,l)*alpha_func(i,j)*e;
+                        for j = 1:obj.dim
+                            alpha{i,j} = d(i,j)*alpha_func(i,k)*e;
                         end
                     end
-
                     o = 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
-                    o = alpha_tot;
             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 points
+        % corresponding to the number of boundary unknowns
         %
         % boundary -- string
         function H = getBoundaryQuadrature(obj, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
 
-            switch boundary
-                case {'w','e'}
-                    j = 1;
-                case {'s','n'}
-                    j = 2;
-            end
-            H = obj.H_boundary{j};
+            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