--- a/+scheme/Elastic2dVariable.m	Sat Jan 26 15:38:58 2019 -0800
+++ b/+scheme/Elastic2dVariable.m	Sat Jan 26 16:56:35 2019 -0800
@@ -549,27 +549,45 @@
         % op -- string
         function o = getBoundaryOperator(obj, op, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
-            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha'})
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'alpha', 'alpha1', 'alpha2'})
 
             switch op
 
                 case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
                     o = obj.([op, '_', boundary]);
 
+                % Yields vector-valued penalty strength given displacement BC on all components
                 case 'alpha'
-                    % alpha = alpha(i,j) is the penalty strength for displacement BC.
-                    e     = obj.getBoundaryOperatorForScalarField('e', boundary);
-                    e_tot = obj.getBoundaryOperator('e', boundary);
-                    alpha = obj.getBoundaryOperatorForScalarField('alpha', boundary);
+                    e               = obj.getBoundaryOperator('e', boundary);
+                    e_scalar        = obj.getBoundaryOperatorForScalarField('e', boundary);
+                    alpha_scalar    = obj.getBoundaryOperatorForScalarField('alpha', boundary);
                     E = obj.E;
-                    [m, n] = size(alpha{1,1});
-                    alpha_tot = sparse(m*obj.dim, n*obj.dim);
+                    [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_tot = alpha_tot + (e_tot'*E{i}*e*alpha{i,l}'*E{l}')';
+                            alpha = alpha + (e'*E{i}*e_scalar*alpha_scalar{i,l}'*E{l}')';
                         end
                     end
-                    o = alpha_tot;
+                    o = alpha;
+
+                % 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);
+
+                    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
 
         end
@@ -586,7 +604,9 @@
                     o = obj.(['e_scalar', '_', boundary]);
 
                 case 'alpha'
-                    % alpha = alpha(i,j) is the penalty strength for displacement BC.
+
+                    % 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;
@@ -610,23 +630,24 @@
 
                     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 j = 1:obj.dim
                             alpha{i,j} = d(i,j)*alpha_func(i,k)*e;
                         end
                     end
-
                     o = alpha;
             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'})