--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Mon Oct 14 13:49:50 2019 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Mon Oct 14 13:54:50 2019 -0700
@@ -24,6 +24,7 @@
         D  % Total operator
 
         Dx, Dy % Physical derivatives
+        n_w, n_e, n_s, n_n % Physical normals
 
         % Boundary operators in cell format, used for BC
         T_w, T_e, T_s, T_n
@@ -47,6 +48,7 @@
         e1_w, e1_e, e1_s, e1_n  % Act on vector field, return scalar field at boundary
         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
 
         % E{i}^T picks out component i
         E
@@ -231,15 +233,15 @@
             obj.e_scalar_s = refObj.e_scalar_s;
             obj.e_scalar_n = refObj.e_scalar_n;
 
-            e1_w = (obj.e_scalar_w'*obj.E{1}')';
-            e1_e = (obj.e_scalar_e'*obj.E{1}')';
-            e1_s = (obj.e_scalar_s'*obj.E{1}')';
-            e1_n = (obj.e_scalar_n'*obj.E{1}')';
+            e1_w = obj.e1_w;
+            e1_e = obj.e1_e;
+            e1_s = obj.e1_s;
+            e1_n = obj.e1_n;
 
-            e2_w = (obj.e_scalar_w'*obj.E{2}')';
-            e2_e = (obj.e_scalar_e'*obj.E{2}')';
-            e2_s = (obj.e_scalar_s'*obj.E{2}')';
-            e2_n = (obj.e_scalar_n'*obj.E{2}')';
+            e2_w = obj.e2_w;
+            e2_e = obj.e2_e;
+            e2_s = obj.e2_s;
+            e2_n = obj.e2_n;
 
             obj.tau1_w = (spdiag(1./s_w)*refObj.tau1_w')';
             obj.tau1_e = (spdiag(1./s_e)*refObj.tau1_e')';
@@ -256,6 +258,53 @@
             obj.tau_s = (refObj.e_s'*obj.e1_s*obj.tau1_s')' + (refObj.e_s'*obj.e2_s*obj.tau2_s')';
             obj.tau_n = (refObj.e_n'*obj.e1_n*obj.tau1_n')' + (refObj.e_n'*obj.e2_n*obj.tau2_n')';
 
+            % Physical normals
+            e_w = obj.e_scalar_w;
+            e_e = obj.e_scalar_e;
+            e_s = obj.e_scalar_s;
+            e_n = obj.e_scalar_n;
+
+            e_w_vec = obj.e_w;
+            e_e_vec = obj.e_e;
+            e_s_vec = obj.e_s;
+            e_n_vec = obj.e_n;
+
+            nu_w = [-1,0];
+            nu_e = [1,0];
+            nu_s = [0,-1];
+            nu_n = [0,1];
+
+            obj.n_w = cell(2,1);
+            obj.n_e = cell(2,1);
+            obj.n_s = cell(2,1);
+            obj.n_n = cell(2,1);
+
+            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);
+
+            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);
+
+            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);
+
+            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);
+
+            % Operators that extract the normal component
+            obj.en_w = (obj.n_w{1}*obj.e1_w' + obj.n_w{2}*obj.e2_w')';
+            obj.en_e = (obj.n_e{1}*obj.e1_e' + obj.n_e{2}*obj.e2_e')';
+            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')';
+
             % Misc.
             obj.refObj = refObj;
             obj.m = refObj.m;
@@ -312,30 +361,6 @@
             [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
         end
 
-        % Returns the component number that is the tangential/normal component
-        % at the specified boundary
-        function comp = getComponent(obj, comp_str, boundary)
-            assertIsMember(comp_str, {'normal', 'tangential'});
-            assertIsMember(boundary, {'w', 'e', 's', 'n'});
-
-            switch boundary
-            case {'w', 'e'}
-                switch comp_str
-                case 'normal'
-                    comp = 1;
-                case 'tangential'
-                    comp = 2;
-                end
-            case {'s', 'n'}
-                switch comp_str
-                case 'normal'
-                    comp = 2;
-                case 'tangential'
-                    comp = 1;
-                end
-            end
-        end
-
         % Returns h11 for the boundary specified by the string boundary.
         % op -- string
         function h11 = getBorrowing(obj, boundary)
@@ -350,31 +375,20 @@
         end
 
         % Returns the outward unit normal vector for the boundary specified by the string boundary.
-        function nu = getNormal(obj, boundary)
+        % n is a cell of diagonal matrices for each normal component, n{1} = n_1, n{2} = n_2.
+        function n = getNormal(obj, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
 
-            switch boundary
-            case 'w'
-                nu = [-1,0];
-            case 'e'
-                nu = [1,0];
-            case 's'
-                nu = [0,-1];
-            case 'n'
-                nu = [0,1];
-            end
+            n = obj.(['n_' 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'})
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en'})
 
-            switch op
-                case {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2'}
-                    o = obj.([op, '_', boundary]);
-            end
+            o = obj.([op, '_', boundary]);
 
         end