--- a/+scheme/Elastic2dStaggeredAnisotropic.m	Sun Jun 07 11:28:22 2020 -0700
+++ b/+scheme/Elastic2dStaggeredAnisotropic.m	Sun Jun 07 11:30:49 2020 -0700
@@ -103,7 +103,11 @@
                     for j = 1:dim
                         for k = 1:dim
                             for l = 1:dim
-                                C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
+                                if numel(C) == nGrids
+                                    C_mat{a}{i,j,k,l} = spdiag(C{a}{i,j,k,l});
+                                else
+                                    C_mat{a}{i,j,k,l} = spdiag(grid.evalOn(g_s{a}, C{i,j,k,l}));
+                                end
                             end
                         end
                     end
@@ -230,7 +234,7 @@
 
 
             % D1_u2s{a, b}{i} approximates ddi and
-            %takes from u grid number b to s grid number a
+            % takes from u grid number b to s grid number a
             % Some of D1_x2y{a, b} are 0.
             D1_u2s = cell(nGrids, nGrids);
             D1_s2u = cell(nGrids, nGrids);
@@ -543,7 +547,7 @@
             % Preallocate
             [~, col] = size(tau{1}{1});
             closure = sparse(m_tot, m_tot);
-            penalty = cell(nGrids, 1);
+            penalty = cell(1, nGrids);
             for a = 1:nGrids
                 [~, col] = size(e_u{a});
                 penalty{a} = sparse(m_tot, col);
@@ -619,6 +623,8 @@
             otherwise
                 error('No such boundary condition: type = %s',type);
             end
+
+            penalty = cell2mat(penalty);
         end
 
         % type     Struct that specifies the interface coupling.
@@ -651,106 +657,127 @@
             v = neighbour_scheme;
 
             % Operators, u side
-            e_u       = u.getBoundaryOperatorForScalarField('e', boundary);
+            eu_u       = u.getBoundaryOperatorForScalarField('e_u', boundary);
+            es_u       = u.getBoundaryOperatorForScalarField('e_s', boundary);
             tau_u     = u.getBoundaryOperator('tau', boundary);
-            h11_u     = u.getBorrowing(boundary);
             nu_u      = u.getNormal(boundary);
 
-            E_u = u.E;
+            G_u = u.G;
+            U_u = u.U;
             C_u = u.C;
-            m_tot_u = u.grid.N();
+            m_tot_u = u.N;
 
             % Operators, v side
-            e_v       = v.getBoundaryOperatorForScalarField('e', neighbour_boundary);
+            eu_v       = v.getBoundaryOperatorForScalarField('e_u', neighbour_boundary);
+            es_v       = v.getBoundaryOperatorForScalarField('e_s', neighbour_boundary);
             tau_v     = v.getBoundaryOperator('tau', neighbour_boundary);
-            h11_v     = v.getBorrowing(neighbour_boundary);
             nu_v      = v.getNormal(neighbour_boundary);
 
-            E_v = v.E;
+            G_v = v.G;
+            U_v = v.U;
             C_v = v.C;
-            m_tot_v = v.grid.N();
-
-            % Fix {'e', 's'}, {'w', 'n'}, and {'x','x'} couplings
-            flipFlag = false;
-            e_v_flip = e_v;
-            if (strcmp(boundary,'s') && strcmp(neighbour_boundary,'e')) || ...
-               (strcmp(boundary,'e') && strcmp(neighbour_boundary,'s')) || ...
-               (strcmp(boundary,'w') && strcmp(neighbour_boundary,'n')) || ...
-               (strcmp(boundary,'n') && strcmp(neighbour_boundary,'w')) || ...
-               (strcmp(boundary,'s') && strcmp(neighbour_boundary,'s')) || ...
-               (strcmp(boundary,'n') && strcmp(neighbour_boundary,'n')) || ...
-               (strcmp(boundary,'w') && strcmp(neighbour_boundary,'w')) || ...
-               (strcmp(boundary,'e') && strcmp(neighbour_boundary,'e'))
-
-                flipFlag = true;
-                e_v_flip = fliplr(e_v);
-
-                t1 = tau_v(:,1:2:end-1);
-                t2 = tau_v(:,2:2:end);
-
-                t1 = fliplr(t1);
-                t2 = fliplr(t2);
-
-                tau_v(:,1:2:end-1) = t1;
-                tau_v(:,2:2:end) = t2;
-            end
+            m_tot_v = v.N;
 
             % Operators that are only required for own domain
-            Hi      = u.Hi_kron;
-            RHOi    = u.RHOi_kron;
-            e_kron  = u.getBoundaryOperator('e', boundary);
+            % Hi      = u.Hi_kron;
+            % RHOi    = u.RHOi_kron;
+            % e_kron  = u.getBoundaryOperator('e', boundary);
+            H       = u.H_u;
+            RHO     = u.RHO;
             T_u     = u.getBoundaryTractionOperator(boundary);
 
             % Shared operators
             H_gamma         = u.getBoundaryQuadratureForScalarField(boundary);
-            H_gamma_kron    = u.getBoundaryQuadrature(boundary);
+            % H_gamma_kron    = u.getBoundaryQuadrature(boundary);
             dim             = u.dim;
+            nGrids          = obj.nGrids;
 
             % 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);
+            % [~, m_int] = size(H_gamma);
+            closure = sparse(m_tot_u, m_tot_u);
+            penalty = sparse(m_tot_u, m_tot_v);
 
-            % ---- Continuity of displacement ------
-
-            % Y: symmetrizing part of penalty
-            % Z: symmetric part of penalty
-            % X = Y + Z.
+            %---- Grid layout -------
+            % gu1 = xp o yp;
+            % gu2 = xd o yd;
+            % gs1 = xd o yp;
+            % gs2 = xp o yd;
+            %------------------------
 
-            % Loop over components to couple across interface
-            for j = 1:dim
+            switch boundary
+                case {'w', 'e'}
+                    switch neighbour_boundary
+                    case {'w', 'e'}
+                        gridCombos = {{1,1,1}, {2,2,2}};
+                    case {'s', 'n'}
+                        gridCombos = {{1,1,2}, {2,2,1}};
+                    end
+                case {'s', 'n'}
+                    switch neighbour_boundary
+                    case {'s', 'n'}
+                        gridCombos = {{2,1,1}, {1,2,2}};
+                    case {'w', 'e'}
+                        gridCombos = {{2,1,2}, {1,2,1}};
+                    end
+            end
 
-                % Loop over components that penalties end up on
-                for i = 1:dim
-                    Y = 1/2*T_u{j,i}';
-                    Z_u = sparse(m_int, m_int);
-                    Z_v = sparse(m_int, m_int);
-                    for l = 1:dim
-                        for k = 1:dim
-                            Z_u = Z_u + e_u'*nu_u(l)*C_u{l,i,k,j}*nu_u(k)*e_u;
-                            Z_v = Z_v + e_v'*nu_v(l)*C_v{l,i,k,j}*nu_v(k)*e_v;
+            % Symmetrizing part
+            for c = 1:nGrids
+                for m = 1:numel(gridCombos)
+                    gc = gridCombos{m};
+                    a = gc{1};
+                    b = gc{2};
+
+                    for i = 1:dim
+                        for j = 1:dim
+                            Y = 1/2*T_u{a,c}{j,i}';
+                            closure = closure + G_u{c}*U_u{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(eu_u{b}'*U_u{b}{j}'*G_u{b}') ));
+                            penalty = penalty - G_u{c}*U_u{c}{i}*((RHO{c}*H{c})\(Y'*H_gamma{a}*(eu_v{b}'*U_v{b}{j}'*G_v{b}') ));
                         end
                     end
-
-                    if flipFlag
-                        Z_v = rot90(Z_v,2);
-                    end
-
-                    Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
-                    X = Y + Z*e_u';
-                    closure = closure + E_u{i}*X'*H_gamma*e_u'*E_u{j}';
-                    penalty = penalty - E_u{i}*X'*H_gamma*e_v_flip'*E_v{j}';
-
                 end
             end
 
-            % ---- Continuity of traction ------
-            closure = closure - 1/2*e_kron*H_gamma_kron*tau_u';
-            penalty = penalty - 1/2*e_kron*H_gamma_kron*tau_v';
+            % Symmetric part
+            for m = 1:numel(gridCombos)
+                gc = gridCombos{m};
+                a = gc{1};
+                b = gc{2};
+                bv = gc{3};
+
+                h11_u = u.getBorrowing(b, boundary);
+                h11_v = v.getBorrowing(bv, neighbour_boundary);
 
-            % ---- Multiply by inverse of density x quadraure ----
-            closure = RHOi*Hi*closure;
-            penalty = RHOi*Hi*penalty;
+                for i = 1:dim
+                    for j = 1:dim
+                        Z_u = 0*es_u{b}'*es_u{b};
+                        Z_v = 0*es_v{bv}'*es_v{bv};
+                        for l = 1:dim
+                            for k = 1:dim
+                                Z_u = Z_u + es_u{b}'*nu_u(l)*C_u{b}{l,i,k,j}*nu_u(k)*es_u{b};
+                                Z_v = Z_v + es_v{bv}'*nu_v(l)*C_v{bv}{l,i,k,j}*nu_v(k)*es_v{bv};
+                            end
+                        end
+                        Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );
+                        % X = es_u{b}'*Z*es_u{b};
+                        % X = Z;
+                        closure = closure + G_u{a}*U_u{a}{i}*((RHO{a}*H{a})\(eu_u{a}*Z'*H_gamma{b}*(eu_u{a}'*U_u{a}{j}'*G_u{a}' ) ));
+                        penalty = penalty - G_u{a}*U_u{a}{i}*((RHO{a}*H{a})\(eu_u{a}*Z'*H_gamma{b}*(eu_v{a}'*U_v{a}{j}'*G_v{a}' ) ));
+                    end
+                end
+            end
+
+            % Continuity of traction
+            for j = 1:dim
+                for m = 1:numel(gridCombos)
+                    gc = gridCombos{m};
+                    a = gc{1};
+                    b = gc{2};
+                    bv = gc{3};
+                    closure = closure - 1/2*G_u{a}*U_u{a}{j}*(RHO{a}\(H{a}\(eu_u{a}*H_gamma{b}*tau_u{b}{j}')));
+                    penalty = penalty - 1/2*G_u{a}*U_u{a}{j}*(RHO{a}\(H{a}\(eu_u{a}*H_gamma{b}*tau_v{bv}{j}')));
+                end
+            end
 
         end
 
@@ -871,7 +898,7 @@
         end
 
         function N = size(obj)
-            N = obj.dim*prod(obj.m);
+            N = length(obj.D);
         end
     end
 end