--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Thu Jul 02 15:00:21 2020 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Thu Jul 02 20:05:50 2020 -0700
@@ -36,6 +36,7 @@
         tau1_w, tau1_e, tau1_s, tau1_n  % Return scalar field
         tau2_w, tau2_e, tau2_s, tau2_n  % Return scalar field
         tau_n_w, tau_n_e, tau_n_s, tau_n_n % Return scalar field
+        tau_t_w, tau_t_e, tau_t_s, tau_t_n % Return scalar field
 
         % Inner products
         H
@@ -320,16 +321,16 @@
             obj.et_s = (obj.tangent_s{1}*obj.e1_s' + obj.tangent_s{2}*obj.e2_s')';
             obj.et_n = (obj.tangent_n{1}*obj.e1_n' + obj.tangent_n{2}*obj.e2_n')';
 
-            % obj.tau_n_w = (obj.en_w'*obj.e_w*obj.tau_w')';
-            % obj.tau_n_e = (obj.en_e'*obj.e_e*obj.tau_e')';
-            % obj.tau_n_s = (obj.en_s'*obj.e_s*obj.tau_s')';
-            % obj.tau_n_n = (obj.en_n'*obj.e_n*obj.tau_n')';
-
             obj.tau_n_w = (obj.n_w{1}*obj.tau1_w' + obj.n_w{2}*obj.tau2_w')';
             obj.tau_n_e = (obj.n_e{1}*obj.tau1_e' + obj.n_e{2}*obj.tau2_e')';
             obj.tau_n_s = (obj.n_s{1}*obj.tau1_s' + obj.n_s{2}*obj.tau2_s')';
             obj.tau_n_n = (obj.n_n{1}*obj.tau1_n' + obj.n_n{2}*obj.tau2_n')';
 
+            obj.tau_t_w = (obj.tangent_w{1}*obj.tau1_w' + obj.tangent_w{2}*obj.tau2_w')';
+            obj.tau_t_e = (obj.tangent_e{1}*obj.tau1_e' + obj.tangent_e{2}*obj.tau2_e')';
+            obj.tau_t_s = (obj.tangent_s{1}*obj.tau1_s' + obj.tangent_s{2}*obj.tau2_s')';
+            obj.tau_t_n = (obj.tangent_n{1}*obj.tau1_n' + obj.tangent_n{2}*obj.tau2_n')';
+
             % Stress operators
             sigma = cell(dim, dim);
             D1 = {obj.Dx, obj.Dy};
@@ -398,17 +399,17 @@
                 end
                 e1 = obj.getBoundaryOperator('e1', boundary);
 
-                bc = {1, type};
-                [c1, p1] = obj.refObj.boundary_condition(boundary, bc, tuning);
-                bc = {2, type};
-                c2 = obj.refObj.boundary_condition(boundary, bc, tuning);
+                bc1 = {1, type};
+                [c1, p1] = obj.refObj.boundary_condition(boundary, bc1, tuning);
+                bc2 = {2, type};
+                c2 = obj.refObj.boundary_condition(boundary, bc2, tuning);
 
                 switch type
                 case {'F','f','Free','free','traction','Traction','t','T'}
                     closure = en*en'*(c1+c2);
                     penalty = en*e1'*p1;
                 case {'D','d','dirichlet','Dirichlet','displacement','Displacement'}
-                    error('Not implemented');
+                    [closure, penalty] = obj.displacementBCNormalTangential(boundary, bc, tuning);
                 end
 
             end
@@ -423,8 +424,8 @@
             end
         end
 
-        function [closure, penalty] = displacementBCNormalTangential(boundary, bc, tuning)
-            error('not implemented');
+        function [closure, penalty] = displacementBCNormalTangential(obj, boundary, bc, tuning)
+            disp('WARNING: displacementBCNormalTangential is only guaranteed to work for displaacement BC on one component and traction on the other');
             u = obj;
 
             component = bc{1};
@@ -432,60 +433,55 @@
 
             switch component
             case 'n'
-                en_u       = u.getBoundaryOperator('en', boundary);
-                tau_n_u     = u.getBoundaryOperator('tau_n', boundary);
+                en      = u.getBoundaryOperator('en', boundary);
+                tau_n   = u.getBoundaryOperator('tau_n', boundary);
+                N       = u.getNormal(boundary);
             case 't'
-                en_u       = u.getBoundaryOperator('et', boundary);
-                tau_n_u     = u.getBoundaryOperator('tau_t', boundary);
+                en      = u.getBoundaryOperator('et', boundary);
+                tau_n   = u.getBoundaryOperator('tau_t', boundary);
+                N       = u.getTangent(boundary);
             end
 
             % Operators
-            e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
-            h11_u     = u.getBorrowing(boundary);
-            n_u      = u.getNormal(boundary);
+            e       = u.getBoundaryOperatorForScalarField('e', boundary);
+            h11     = u.getBorrowing(boundary);
+            n      = u.getNormal(boundary);
 
-            C_u = u.C;
-            Ji_u = u.Ji;
-            s_u = spdiag(u.(['s_' boundary]));
-            m_tot_u = u.grid.N();
+            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}';
 
-            % Shared operators
             H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
             dim             = u.dim;
 
             % Preallocate
             [~, m_int] = size(H_gamma);
-            closure = sparse(dim*m_tot_u, dim*m_tot_u);
-            penalty = sparse(dim*m_tot_u, dim*m_tot_v);
+            closure = sparse(dim*m_tot, dim*m_tot);
+            penalty = sparse(dim*m_tot, m_int);
 
-            % Continuity of normal displacement, term 1: The symmetric term
-            Z_u = sparse(m_int, m_int);
-            Z_v = sparse(m_int, 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_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{l,i,k,j}*e_u*n_u{k}*n_u{l};
-                            Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{l,i,k,j}*e_v*n_v{k}*n_v{l};
+                            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/(4*h11_u)*s_u*Z_u + 1/(4*h11_v)*s_v*Z_v );
-            closure = closure + en_u*H_gamma*Z*en_u';
-            penalty = penalty + en_u*H_gamma*Z*en_v';
+            Z = -tuning*dim*1/h11*s*Z;
+            closure = closure + en*H_gamma*Z*en';
+            penalty = penalty - en*H_gamma*Z;
 
-            % Continuity of normal displacement, term 2: The symmetrizing term
-            closure = closure + 1/2*tau_n_u*H_gamma*en_u';
-            penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
-
-            % Continuity of normal traction
-            closure = closure - 1/2*en_u*H_gamma*tau_n_u';
-            penalty = penalty - 1/2*en_u*H_gamma*tau_n_v';
+            % 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;
@@ -565,8 +561,8 @@
                 for j = 1:dim
                     for l = 1:dim
                         for k = 1:dim
-                            Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{l,i,k,j}*e_u*n_u{k}*n_u{l};
-                            Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{l,i,k,j}*e_v*n_v{k}*n_v{l};
+                            Z_u = Z_u + n_u{i}*n_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*n_u{l};
+                            Z_v = Z_v + n_v{i}*n_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*n_v{l};
                         end
                     end
                 end
@@ -615,11 +611,19 @@
             n = obj.(['n_' boundary]);
         end
 
+        % Returns the unit tangent vector for the boundary specified by the string boundary.
+        % t is a cell of diagonal matrices for each normal component, t{1} = t_1, t{2} = t_2.
+        function t = getTangent(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
+            t = obj.(['tangent_' boundary]);
+        end
+
         % Returns the boundary operator op for the boundary specified by the string boundary.
         % op -- string
         function o = getBoundaryOperator(obj, op, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
-            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n'})
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n', 'tau_t'})
 
             o = obj.([op, '_', boundary]);