--- a/+scheme/Elastic2dVariableAnisotropic.m	Thu Jul 02 20:05:50 2020 -0700
+++ b/+scheme/Elastic2dVariableAnisotropic.m	Fri Jul 10 20:24:23 2020 -0700
@@ -55,7 +55,8 @@
 
         % The coefficients can either be function handles or grid functions
         % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
-        function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag)
+        function obj = Elastic2dVariableAnisotropic(g, order, rho, C, opSet, optFlag, hollow)
+            default_arg('hollow', false);
             default_arg('rho', @(x,y) 0*x+1);
             default_arg('opSet',{@sbp.D2VariableCompatible, @sbp.D2VariableCompatible});
             default_arg('optFlag', false);
@@ -196,39 +197,66 @@
             switch order
             case 2
                 width = 3;
+                nBP = 2;
             case 4
                 width = 5;
+                nBP = 6;
             case 6
                 width = 7;
+                nBP = 9;
             end
             for j = 1:dim
                 for k = 1:dim
                     for l = 1:dim
-                        D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
+                        if hollow
+                            D2_temp{j,k,l} = sparse(m_tot, m_tot);
+                        else
+                            D2_temp{j,k,l} = spalloc(m_tot, m_tot, width*m_tot);
+                        end
                     end
                 end
             end
             ind = grid.funcToMatrix(g, 1:m_tot);
 
             k = 1;
+            if hollow
+                mask = sparse(m(1), m(1));
+                mask(1:nBP, 1:nBP) = speye(nBP, nBP);
+                mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP);
+                maskX = kron(mask, speye(m(2), m(2)));
+                maskX = E{1}*maskX*E{1}' + E{2}*maskX*E{2}';
+            end
             for r = 1:m(2)
                 p = ind(:,r);
                 for j = 1:dim
                     for l = 1:dim
                         coeff = C{k,j,k,l};
                         D_kk = D2{1}(coeff(p));
+                        if hollow && r > nBP && r < m(2) - nBP + 1
+                            D_kk = mask*D_kk;
+                        end
                         D2_temp{j,k,l}(p,p) = D_kk;
                     end
                 end
             end
 
             k = 2;
+            if hollow
+                mask = sparse(m(2), m(2));
+                mask(1:nBP, 1:nBP) = speye(nBP, nBP);
+                mask(end-nBP+1:end, end-nBP+1:end) = speye(nBP, nBP);
+                maskY = kron(speye(m(1), m(1)), mask);
+                maskY = E{1}*maskY*E{1}' + E{2}*maskY*E{2}';
+            end
             for r = 1:m(1)
                 p = ind(r,:);
                 for j = 1:dim
                     for l = 1:dim
                         coeff = C{k,j,k,l};
                         D_kk = D2{2}(coeff(p));
+                        if hollow && r > nBP && r < m(1) - nBP + 1
+                            D_kk = mask*D_kk;
+                        end
                         D2_temp{j,k,l}(p,p) = D_kk;
                     end
                 end
@@ -261,6 +289,13 @@
                 end
             end
             clear D2_temp;
+            if hollow
+                mask = maskX + maskY;
+                mask = mask>0;
+
+                % D = maskX*maskY*D;
+                D = mask*D;
+            end
             D = obj.RHOi_kron*D;
             obj.D = D;
             clear D;