--- a/+scheme/ViscoElastic2d.m	Mon Jul 20 17:45:58 2020 -0700
+++ b/+scheme/ViscoElastic2d.m	Mon Jul 20 22:06:33 2020 -0700
@@ -292,11 +292,13 @@
                 end
 
                 % Add viscous part of closure
-                for j = 1:dim
-                    for i = 1:dim
+                for i = 1:dim
+                    for j = 1:dim
                         for k = 1:dim
                             for l = 1:dim
-                                closure = closure + eU{j}*( (RHO*H)\(C{i,j,k,l}*e*dir{j}^2*H_gamma*n{i}*e'*eGamma{k,l}') );
+                                for m = 1:dim
+                                    closure = closure + eU{m}*( (RHO*H)\(C{i,j,k,l}*e*dir{m}*dir{j}*H_gamma*n{i}*e'*eGamma{k,l}') );
+                                end
                             end
                         end
                     end
@@ -309,13 +311,27 @@
                 closure = Eu*closure*Eu';
                 penalty = Eu*penalty;
 
+                switch component
+                case 't'
+                    dir = obj.getTangent(boundary);
+                case 'n'
+                    dir = obj.getNormal(boundary);
+                case 1
+                    dir = {1, 0};
+                case 2
+                    dir = {0, 1};
+                end
+
                 % Add penalty to strain rate eq
-                l = component;
                 for i = 1:dim
                     for j = 1:dim
                         for k = 1:dim
-                            closure = closure - eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}*e'*eU{l}') );
-                            penalty = penalty + eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*n{k}) );
+                            for l = 1:dim
+                                for m = 1:dim
+                                    closure = closure - eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*dir{l}*dir{m}*n{k}*e'*eU{m}') );
+                                end
+                                penalty = penalty + eGamma{i,j}*( (H*ETA)\(C{i,j,k,l}*e*H_gamma*dir{l}*n{k}) );
+                            end
                         end
                     end
                 end
@@ -324,70 +340,6 @@
 
         end
 
-        function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
-            disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displacement BC on one component and traction on the other');
-            u = obj;
-
-            component = bc{1};
-            type = bc{2};
-
-            switch component
-            case 'n'
-                en      = u.getBoundaryOperator('en', boundary);
-                tau_n   = u.getBoundaryOperator('tau_n', boundary);
-                N       = u.getNormal(boundary);
-            case 't'
-                en      = u.getBoundaryOperator('et', boundary);
-                tau_n   = u.getBoundaryOperator('tau_t', boundary);
-                N       = u.getTangent(boundary);
-            end
-
-            % Operators
-            e       = u.getBoundaryOperatorForScalarField('e', boundary);
-            h11     = u.getBorrowing(boundary);
-            n      = u.getNormal(boundary);
-
-            C = u.C;
-            Ji = u.Ji;
-            s = spdiag(u.(['s_' boundary]));
-            m_tot = u.grid.N();
-
-            Hi      = u.E{1}*inv(u.H)*u.E{1}' + u.E{2}*inv(u.H)*u.E{2}';
-            RHOi    = u.E{1}*inv(u.RHO)*u.E{1}' + u.E{2}*inv(u.RHO)*u.E{2}';
-
-            H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
-            dim             = u.dim;
-
-            % Preallocate
-            [~, m_int] = size(H_gamma);
-            closure = sparse(dim*m_tot, dim*m_tot);
-            penalty = sparse(dim*m_tot, m_int);
-
-            % Term 1: The symmetric term
-            Z = sparse(m_int, m_int);
-            for i = 1:dim
-                for j = 1:dim
-                    for l = 1:dim
-                        for k = 1:dim
-                            Z = Z + n{i}*N{j}*e'*Ji*C{i,j,k,l}*e*n{k}*N{l};
-                        end
-                    end
-                end
-            end
-
-            Z = -tuning*dim*1/h11*s*Z;
-            closure = closure + en*H_gamma*Z*en';
-            penalty = penalty - en*H_gamma*Z;
-
-            % Term 2: The symmetrizing term
-            closure = closure + tau_n*H_gamma*en';
-            penalty = penalty - tau_n*H_gamma;
-
-            % Multiply all normal component terms by inverse of density x quadraure
-            closure = RHOi*Hi*closure;
-            penalty = RHOi*Hi*penalty;
-        end
-
         % type     Struct that specifies the interface coupling.
         %          Fields:
         %          -- tuning:           penalty strength, defaults to 1.0