--- a/+scheme/Elastic2dVariableAnisotropic.m	Sat Sep 07 13:05:26 2019 -0700
+++ b/+scheme/Elastic2dVariableAnisotropic.m	Fri Sep 20 14:57:02 2019 -0700
@@ -140,9 +140,11 @@
             end
 
             %====== Assemble full operators ========
+            I_dim = speye(dim, dim);
             RHO = spdiag(rho);
             obj.RHO = RHO;
             obj.RHOi = inv(RHO);
+            obj.RHOi_kron = kron(obj.RHOi, I_dim);
 
             obj.D1 = cell(dim,1);
             obj.D2 = cell(dim,dim,dim);
@@ -164,7 +166,6 @@
             e_scalar_s = e_l{2};
             e_scalar_n = e_r{2};
 
-            I_dim = speye(dim, dim);
             e_w = kron(e_scalar_w, I_dim);
             e_e = kron(e_scalar_e, I_dim);
             e_s = kron(e_scalar_s, I_dim);
@@ -192,10 +193,10 @@
 
 
             % D2
-            for i = 1:dim
-                for k = 2:dim
-                    for l = 2:dim
-                        obj.D2{i,k,l} = sparse(m_tot);
+            for j = 1:dim
+                for k = 1:dim
+                    for l = 1:dim
+                        obj.D2{j,k,l} = sparse(m_tot,m_tot);
                     end
                 end
             end
@@ -204,11 +205,11 @@
             k = 1;
             for r = 1:m(2)
                 p = ind(:,r);
-                for i = 1:dim
+                for j = 1:dim
                     for l = 1:dim
-                        coeff = C{i,k,k,l};
+                        coeff = C{k,j,k,l};
                         D_kk = D2{1}(coeff(p));
-                        obj.D2{i,k,l}(p,p) = D_kk;
+                        obj.D2{j,k,l}(p,p) = D_kk;
                     end
                 end
             end
@@ -216,11 +217,11 @@
             k = 2;
             for r = 1:m(1)
                 p = ind(r,:);
-                for i = 1:dim
+                for j = 1:dim
                     for l = 1:dim
-                        coeff = C{i,k,k,l};
+                        coeff = C{k,j,k,l};
                         D_kk = D2{2}(coeff(p));
-                        obj.D2{i,k,l}(p,p) = D_kk;
+                        obj.D2{j,k,l}(p,p) = D_kk;
                     end
                 end
             end
@@ -238,19 +239,20 @@
             D2 = obj.D2;
             D1 = obj.D1;
             D = sparse(dim*m_tot,dim*m_tot);
-            d = @kroneckerDelta;    % Kronecker delta
-            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
             for i = 1:dim
                 for j = 1:dim
                     for k = 1:dim
                         for l = 1:dim
-                            D = D + E{i}*inv(RHO)*( d(j,k)*D2{i,k,l}*E{l}' +...
-                                                    db(j,k)*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}' ...
-                                                  );
+                            if i == k
+                                D = D + E{j}*D2{j,k,l}*E{l}';
+                            else
+                                D = D + E{j}*D1{i}*C_mat{i,j,k,l}*D1{k}*E{l}';
+                            end
                         end
                     end
                 end
             end
+            D = obj.RHOi_kron*D;
             obj.D = D;
             %=========================================%'
 
@@ -288,9 +290,9 @@
                         % Derivative direction k
                         for k = 1:dim
                             T_l{j}{i,l} = T_l{j}{i,l} ...
-                                        - (e_l{j}'*C_mat{i,j,k,l}*D1{k})';
+                                        - (e_l{j}'*C_mat{j,i,k,l}*D1{k})';
                             T_r{j}{i,l} = T_r{j}{i,l} ...
-                                        + (e_r{j}'*C_mat{i,j,k,l}*D1{k})';
+                                        + (e_r{j}'*C_mat{j,i,k,l}*D1{k})';
                         end
                         tau_l{j}{i} = tau_l{j}{i} + (T_l{j}{i,l}'*E{l}')';
                         tau_r{j}{i} = tau_r{j}{i} + (T_r{j}{i,l}'*E{l}')';
@@ -344,7 +346,6 @@
             obj.tau_n = (e_n'*e1_n*obj.tau1_n')' + (e_n'*e2_n*obj.tau2_n')';
 
             % Kroneckered norms and coefficients
-            obj.RHOi_kron = kron(obj.RHOi, I_dim);
             obj.Hi_kron = kron(obj.Hi, I_dim);
 
             % Misc.
@@ -435,7 +436,7 @@
                     Z = sparse(m_tot, m_tot);
                     for l = 1:dim
                         for k = 1:dim
-                            Z = Z + nu(l)*C{i,l,k,j}*nu(k);
+                            Z = Z + nu(l)*C{l,i,k,j}*nu(k);
                         end
                     end
                     Z = -tuning*dim/h11*Z;
@@ -536,8 +537,8 @@
                     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{i,l,k,j}*nu_u(k)*e_u;
-                            Z_v = Z_v + e_v'*nu_v(l)*C_v{i,l,k,j}*nu_v(k)*e_v;
+                            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;
                         end
                     end
                     Z = -tuning*dim*( 1/(4*h11_u)*Z_u + 1/(4*h11_v)*Z_v );