--- a/+scheme/Elastic2dVariable.m	Fri Sep 14 10:41:00 2018 -0700
+++ b/+scheme/Elastic2dVariable.m	Sun Sep 16 18:00:14 2018 -0700
@@ -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,14 +30,14 @@
         T_l, T_r
         tau_l, tau_r
 
-        H, Hi % Inner products
+        H, Hi, H_1D % Inner products
         phi % Borrowing constant for (d1 - e^T*D1) from R
         gamma % Borrowing constant for d1 from M
         H11 % First element of H
         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
@@ -171,6 +171,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);
@@ -230,14 +231,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}';
@@ -270,7 +271,7 @@
         % 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.
         %       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.
@@ -317,20 +318,20 @@
                 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 ); 
+                                      + 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*H_gamma*(e'*E{k}' ); 
+                    closure = closure + E{i}*RHOi*Hi*B'*e*H_gamma*(e'*E{k}' );
                     penalty = penalty - E{i}*RHOi*Hi*B'*e*H_gamma;
-                end 
+                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* (e'*tau{k} );
                     penalty = penalty + E{k}*RHOi*Hi*e*H_gamma;
 
             % Unknown boundary condition
@@ -357,7 +358,7 @@
             Hi = obj.Hi;
             RHOi = obj.RHOi;
             dim = obj.dim;
-        
+
             %--- Other operators ----
             m_tot_u = obj.grid.N();
             E = obj.E;
@@ -373,7 +374,7 @@
             lambda_v = e_v'*LAMBDA_v*e_v;
             mu_v = e_v'*MU_v*e_v;
             %-------------------------
-            
+
             % Borrowing constants
             phi_u = obj.phi{j};
             h_u = obj.h(j);
@@ -386,7 +387,7 @@
             gamma_v = neighbour_scheme.gamma{j_v};
 
             % E > sum_i 1/(2*alpha_ij)*(tau_i)^2
-            function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu) 
+            function [alpha_ii, alpha_ij] = computeAlpha(phi,h,h11,gamma,lambda,mu)
                 th1 = h11/(2*dim);
                 th2 = h11*phi/2;
                 th3 = h*gamma;
@@ -398,7 +399,7 @@
             end
 
             [alpha_ii_u, alpha_ij_u] = computeAlpha(phi_u,h_u,h11_u,gamma_u,lambda_u,mu_u);
-            [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);  
+            [alpha_ii_v, alpha_ij_v] = computeAlpha(phi_v,h_v,h11_v,gamma_v,lambda_v,mu_v);
             sigma_ii = tuning*(alpha_ii_u + alpha_ii_v)/4;
             sigma_ij = tuning*(alpha_ij_u + alpha_ij_v)/4;
 
@@ -420,9 +421,9 @@
 
                 % 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 
+                    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
 
@@ -494,7 +495,7 @@
                                 varargout{i} = obj.tau_l{j};
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 varargout{i} = obj.tau_r{j};
-                        end                        
+                        end
                     otherwise
                         error(['No such operator: operator = ' op{i}]);
                 end