--- a/+scheme/Elastic2dVariableAnisotropic.m	Sat Sep 07 13:05:26 2019 -0700
+++ b/+scheme/Elastic2dVariableAnisotropic.m	Sun Sep 15 14:39:42 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);
@@ -193,9 +194,9 @@
 
             % D2
             for i = 1:dim
-                for k = 2:dim
-                    for l = 2:dim
-                        obj.D2{i,k,l} = sparse(m_tot);
+                for k = 1:dim
+                    for l = 1:dim
+                        obj.D2{i,k,l} = sparse(m_tot,m_tot);
                     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 j == k
+                                D = D + E{i}*D2{i,k,l}*E{l}';
+                            else
+                                D = D + E{i}*D1{j}*C_mat{i,j,k,l}*D1{k}*E{l}';
+                            end
                         end
                     end
                 end
             end
+            D = obj.RHOi_kron*D;
             obj.D = D;
             %=========================================%'
 
@@ -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.