--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Sat Jul 25 15:56:31 2020 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Sun Jul 26 17:38:28 2020 -0700
@@ -486,7 +486,7 @@
             closure = closure + tau_n*H_gamma*en';
             penalty = penalty - tau_n*H_gamma;
 
-            % Multiply all normal component terms by inverse of density x quadraure
+            % Multiply all terms by inverse of density x quadraure
             closure = RHOi*Hi*closure;
             penalty = RHOi*Hi*penalty;
         end
@@ -505,6 +505,8 @@
             switch type.type
             case 'standard'
                 [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
+            case 'normalTangential'
+                [closure, penalty] = obj.interfaceNormalTangential(boundary,neighbour_scheme,neighbour_boundary,type);
             case 'frictionalFault'
                 [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
             end
@@ -592,6 +594,122 @@
 
         end
 
+        % Same interface conditions as in interfaceStandard, but imposed in the normal-tangential
+        % coordinate system. For the isotropic case, the components decouple nicely.
+        % The resulting scheme is not identical to that of interfaceStandard. This appears to be better.
+        function [closure, penalty] = interfaceNormalTangential(obj,boundary,neighbour_scheme,neighbour_boundary,type)
+            tuning = type.tuning;
+
+            % u denotes the solution in the own domain
+            % v denotes the solution in the neighbour domain
+
+            u = obj;
+            v = neighbour_scheme;
+
+            % Operators, u side
+            e_u         = u.getBoundaryOperatorForScalarField('e', boundary);
+            en_u        = u.getBoundaryOperator('en', boundary);
+            et_u        = u.getBoundaryOperator('et', boundary);
+            tau_n_u     = u.getBoundaryOperator('tau_n', boundary);
+            tau_t_u     = u.getBoundaryOperator('tau_t', boundary);
+            h11_u       = u.getBorrowing(boundary);
+            n_u         = u.getNormal(boundary);
+            t_u         = u.getTangent(boundary);
+
+            C_u = u.C;
+            Ji_u = u.Ji;
+            s_u = spdiag(u.(['s_' boundary]));
+            m_tot_u = u.grid.N();
+
+            % Operators, v side
+            e_v         = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
+            en_v        = v.getBoundaryOperator('en', neighbour_boundary);
+            et_v        = v.getBoundaryOperator('et', neighbour_boundary);
+            tau_n_v     = v.getBoundaryOperator('tau_n', neighbour_boundary);
+            tau_t_v     = v.getBoundaryOperator('tau_t', neighbour_boundary);
+            h11_v       = v.getBorrowing(neighbour_boundary);
+            n_v         = v.getNormal(neighbour_boundary);
+            t_v         = v.getTangent(neighbour_boundary);
+
+            C_v = v.C;
+            Ji_v = v.Ji;
+            s_v = spdiag(v.(['s_' neighbour_boundary]));
+            m_tot_v = v.grid.N();
+
+            % Operators that are only required for own domain
+            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);
+
+            % -- Continuity of displacement, term 1: The symmetric term ---
+            Zn_u = sparse(m_int, m_int);
+            Zn_v = sparse(m_int, m_int);
+            Zt_u = sparse(m_int, m_int);
+            Zt_v = sparse(m_int, m_int);
+            for i = 1:dim
+                for j = 1:dim
+                    for l = 1:dim
+                        for k = 1:dim
+                            % Penalty strength for normal component
+                            Zn_u = Zn_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};
+                            Zn_v = Zn_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};
+
+                            % Penalty strength for tangential component
+                            Zt_u = Zt_u + n_u{i}*t_u{j}*e_u'*Ji_u*C_u{i,j,k,l}*e_u*n_u{k}*t_u{l};
+                            Zt_v = Zt_v + n_v{i}*t_v{j}*e_v'*Ji_v*C_v{i,j,k,l}*e_v*n_v{k}*t_v{l};
+                        end
+                    end
+                end
+            end
+
+            Zn = -tuning*dim*( 1/(4*h11_u)*s_u*Zn_u + 1/(4*h11_v)*s_v*Zn_v );
+            Zt = -tuning*dim*( 1/(4*h11_u)*s_u*Zt_u + 1/(4*h11_v)*s_v*Zt_v );
+
+            % Continuity of normal component
+            closure = closure + en_u*H_gamma*Zn*en_u';
+            penalty = penalty + en_u*H_gamma*Zn*en_v';
+
+            % Continuity of tangential component
+            closure = closure + et_u*H_gamma*Zt*et_u';
+            penalty = penalty + et_u*H_gamma*Zt*et_v';
+            %------------------------------------------------------------------
+
+            % --- Continuity of displacement, term 2: The symmetrizing term
+
+            % Continuity of normal displacement
+            closure = closure + 1/2*tau_n_u*H_gamma*en_u';
+            penalty = penalty + 1/2*tau_n_u*H_gamma*en_v';
+
+            % Continuity of tangential displacement
+            closure = closure + 1/2*tau_t_u*H_gamma*et_u';
+            penalty = penalty + 1/2*tau_t_u*H_gamma*et_v';
+            % ------------------------------------------------------------------
+
+            % --- Continuity of tractions -----------------------------
+
+            % 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';
+
+            % Continuity of tangential traction
+            closure = closure - 1/2*et_u*H_gamma*tau_t_u';
+            penalty = penalty + 1/2*et_u*H_gamma*tau_t_v';
+            %--------------------------------------------------------------------
+
+            % Multiply all terms by inverse of density x quadraure
+            closure = RHOi*Hi*closure;
+            penalty = RHOi*Hi*penalty;
+
+        end
+
 
         % Returns h11 for the boundary specified by the string boundary.
         % op -- string