--- a/+scheme/Elastic2dVariable.m	Tue Dec 25 07:21:19 2018 +0100
+++ b/+scheme/Elastic2dVariable.m	Tue Dec 25 07:23:38 2018 +0100
@@ -367,7 +367,7 @@
 
             % 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;
@@ -389,7 +389,7 @@
             % Dirichlet boundary condition
             case {'D','d','dirichlet','Dirichlet'}
 
-                alpha = obj.get_boundary_operator('alpha', boundary);
+                alpha = obj.getBoundaryOperator('alpha', boundary);
 
                 % Loop over components that Dirichlet penalties end up on
                 for i = 1:dim
@@ -443,8 +443,8 @@
             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, T, tau, H_gamma] = obj.getBoundaryOperator({'e','T','tau','H'}, boundary);
+            [e_v, tau_v] = neighbour_scheme.getBoundaryOperator({'e','tau'}, neighbour_boundary);
 
             % Operators and quantities that correspond to the own domain only
             Hi = obj.Hi;
@@ -536,7 +536,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'}
@@ -560,6 +561,12 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 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'}
@@ -567,6 +574,7 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 varargout{k} = obj.d1_r{j};
                         end
+
                     case 'T'
                         switch boundary
                             case {'w','W','west','West','s','S','south','South'}
@@ -574,6 +582,7 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 varargout{k} = obj.T_r{j};
                         end
+
                     case 'tau'
                         switch boundary
                             case {'w','W','west','West','s','S','south','South'}
@@ -581,16 +590,25 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 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.
-                        switch boundary
-                            case {'w','W','west','West','s','S','south','South'}
-                                e = obj.e_l{j};
-                            case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
-                                e = obj.e_r{j};
-                        end
+                        e = obj.getBoundaryOperator('e', boundary);
 
                         tuning = 1.2;
                         LAMBDA = obj.LAMBDA;
@@ -619,6 +637,20 @@
                         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{k}]);
                 end
@@ -626,6 +658,20 @@
 
         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