--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Sun Jun 07 11:33:44 2020 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Mon Jun 29 15:23:01 2020 -0700
@@ -26,6 +26,7 @@
         Dx, Dy % Physical derivatives
         sigma % Cell matrix of physical stress operators
         n_w, n_e, n_s, n_n % Physical normals
+        tangent_w, tangent_e, tangent_s, tangent_n % Physical tangents
 
         % Boundary operators in cell format, used for BC
         T_w, T_e, T_s, T_n
@@ -50,6 +51,7 @@
         e2_w, e2_e, e2_s, e2_n  % Act on vector field, return scalar field at boundary
         e_scalar_w, e_scalar_e, e_scalar_s, e_scalar_n; % Act on scalar field, return scalar field
         en_w, en_e, en_s, en_n  % Act on vector field, return normal component
+        et_w, et_e, et_s, et_n  % Act on vector field, return tangential component
 
         % E{i}^T picks out component i
         E
@@ -280,25 +282,30 @@
             obj.n_s = cell(2,1);
             obj.n_n = cell(2,1);
 
+            % Compute normal and rotate (exactly!) 90 degrees counter-clockwise to get tangent
             n_w_1 = (1./s_w).*e_w'*(J.*(K{1,1}*nu_w(1) + K{1,2}*nu_w(2)));
             n_w_2 = (1./s_w).*e_w'*(J.*(K{2,1}*nu_w(1) + K{2,2}*nu_w(2)));
             obj.n_w{1} = spdiag(n_w_1);
             obj.n_w{2} = spdiag(n_w_2);
+            obj.tangent_w = {-obj.n_w{2}, obj.n_w{1}};
 
             n_e_1 = (1./s_e).*e_e'*(J.*(K{1,1}*nu_e(1) + K{1,2}*nu_e(2)));
             n_e_2 = (1./s_e).*e_e'*(J.*(K{2,1}*nu_e(1) + K{2,2}*nu_e(2)));
             obj.n_e{1} = spdiag(n_e_1);
             obj.n_e{2} = spdiag(n_e_2);
+            obj.tangent_e = {-obj.n_e{2}, obj.n_e{1}};
 
             n_s_1 = (1./s_s).*e_s'*(J.*(K{1,1}*nu_s(1) + K{1,2}*nu_s(2)));
             n_s_2 = (1./s_s).*e_s'*(J.*(K{2,1}*nu_s(1) + K{2,2}*nu_s(2)));
             obj.n_s{1} = spdiag(n_s_1);
             obj.n_s{2} = spdiag(n_s_2);
+            obj.tangent_s = {-obj.n_s{2}, obj.n_s{1}};
 
             n_n_1 = (1./s_n).*e_n'*(J.*(K{1,1}*nu_n(1) + K{1,2}*nu_n(2)));
             n_n_2 = (1./s_n).*e_n'*(J.*(K{2,1}*nu_n(1) + K{2,2}*nu_n(2)));
             obj.n_n{1} = spdiag(n_n_1);
             obj.n_n{2} = spdiag(n_n_2);
+            obj.tangent_n = {-obj.n_n{2}, obj.n_n{1}};
 
             % Operators that extract the normal component
             obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
@@ -306,6 +313,12 @@
             obj.en_s = (obj.n_s{1}*obj.e1_s' + obj.n_s{2}*obj.e2_s')';
             obj.en_n = (obj.n_n{1}*obj.e1_n' + obj.n_n{2}*obj.e2_n')';
 
+            % Operators that extract the tangential component
+            obj.et_w = (obj.tangent_w{1}*obj.e1_w' + obj.tangent_w{2}*obj.e2_w')';
+            obj.et_e = (obj.tangent_e{1}*obj.e1_e' + obj.tangent_e{2}*obj.e2_e')';
+            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')';
+
             % Stress operators
             sigma = cell(dim, dim);
             D1 = {obj.Dx, obj.Dy};
@@ -354,15 +367,48 @@
             default_arg('tuning', 1.0);
             assert( iscell(bc), 'The BC type must be a 2x1 or 3x1 cell array' );
 
-            [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
+            component = bc{1};
+            type = bc{2};
+
+            switch component
+
+            % If conditions on Cartesian components
+            case {1,2}
+                [closure, penalty] = obj.refObj.boundary_condition(boundary, bc, tuning);
+
+            % If conditions on normal or tangential components
+            case {'n', 't'}
 
-            type = bc{2};
+                switch component
+                    case 'n'
+                        en = obj.getBoundaryOperator('en', boundary);
+                    case 't'
+                        en = obj.getBoundaryOperator('et', boundary);
+                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);
+
+                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');
+                end
+
+            end
 
             switch type
             case {'F','f','Free','free','traction','Traction','t','T'}
+
                 s = obj.(['s_' boundary]);
                 s = spdiag(s);
                 penalty = penalty*s;
+
             end
         end
 
@@ -404,7 +450,7 @@
         % op -- string
         function o = getBoundaryOperator(obj, op, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
-            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et'})
 
             o = obj.([op, '_', boundary]);