--- a/+scheme/Elastic2dCurvilinearAnisotropic.m	Mon Jun 29 15:26:36 2020 -0700
+++ b/+scheme/Elastic2dCurvilinearAnisotropic.m	Thu Jul 02 15:00:21 2020 -0700
@@ -35,6 +35,7 @@
         tau_w, tau_e, tau_s, tau_n      % Return vector field
         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
 
         % Inner products
         H
@@ -319,6 +320,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')';
+
             % Stress operators
             sigma = cell(dim, dim);
             D1 = {obj.Dx, obj.Dy};
@@ -412,6 +423,75 @@
             end
         end
 
+        function [closure, penalty] = displacementBCNormalTangential(boundary, bc, tuning)
+            error('not implemented');
+            u = obj;
+
+            component = bc{1};
+            type = bc{2};
+
+            switch component
+            case 'n'
+                en_u       = u.getBoundaryOperator('en', boundary);
+                tau_n_u     = u.getBoundaryOperator('tau_n', boundary);
+            case 't'
+                en_u       = u.getBoundaryOperator('et', boundary);
+                tau_n_u     = u.getBoundaryOperator('tau_t', boundary);
+            end
+
+            % Operators
+            e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
+            h11_u     = u.getBorrowing(boundary);
+            n_u      = u.getNormal(boundary);
+
+            C_u = u.C;
+            Ji_u = u.Ji;
+            s_u = spdiag(u.(['s_' boundary]));
+            m_tot_u = 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);
+
+            % Continuity of normal displacement, term 1: The symmetric term
+            Z_u = sparse(m_int, m_int);
+            Z_v = 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};
+                        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';
+
+            % 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';
+
+            % 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
@@ -420,11 +500,100 @@
 
             defaultType.tuning = 1.0;
             defaultType.interpolation = 'none';
+            defaultType.type = 'standard';
             default_struct('type', defaultType);
 
-            [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
+            switch type.type
+            case 'standard'
+                [closure, penalty] = obj.refObj.interface(boundary,neighbour_scheme.refObj,neighbour_boundary,type);
+            case 'frictionalFault'
+                [closure, penalty] = obj.interfaceFrictionalFault(boundary,neighbour_scheme,neighbour_boundary,type);
+            end
+
         end
 
+        function [closure, penalty] = interfaceFrictionalFault(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);
+            tau_n_u     = u.getBoundaryOperator('tau_n', boundary);
+            h11_u     = u.getBorrowing(boundary);
+            n_u      = u.getNormal(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);
+            tau_n_v     = v.getBoundaryOperator('tau_n', neighbour_boundary);
+            h11_v     = v.getBorrowing(neighbour_boundary);
+            n_v      = v.getNormal(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 normal displacement, term 1: The symmetric term
+            Z_u = sparse(m_int, m_int);
+            Z_v = 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};
+                        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';
+
+            % 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';
+
+            % Multiply all normal component terms by inverse of density x quadraure
+            closure = RHOi*Hi*closure;
+            penalty = RHOi*Hi*penalty;
+
+            % ---- Tangential tractions are imposed just like traction BC ------
+            closure = closure + obj.boundary_condition(boundary, {'t', 't'});
+
+        end
+
+
         % Returns h11 for the boundary specified by the string boundary.
         % op -- string
         function h11 = getBorrowing(obj, boundary)
@@ -450,7 +619,7 @@
         % op -- string
         function o = getBoundaryOperator(obj, op, boundary)
             assertIsMember(boundary, {'w', 'e', 's', 'n'})
-            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et'})
+            assertIsMember(op, {'e', 'e1', 'e2', 'tau', 'tau1', 'tau2', 'en', 'et', 'tau_n'})
 
             o = obj.([op, '_', boundary]);