--- a/+scheme/Elastic2dVariable.m	Sat May 04 14:44:34 2019 -0700
+++ b/+scheme/Elastic2dVariable.m	Sun May 05 19:05:31 2019 -0700
@@ -61,11 +61,13 @@
     methods
 
         % The coefficients can either be function handles or grid functions
-        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
+        % optFlag -- if true, extra computations are performed, which may be helpful for optimization.
+        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet, optFlag)
             default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
             default_arg('lambda', @(x,y) 0*x+1);
             default_arg('mu', @(x,y) 0*x+1);
             default_arg('rho', @(x,y) 0*x+1);
+            default_arg('optFlag', false);
             dim = 2;
 
             assert(isa(g, 'grid.Cartesian'))
@@ -356,42 +358,46 @@
             obj.dim = dim;
 
             % B, used for adjoint optimization
-            B = cell(dim, 1);
-            for i = 1:dim
-                B{i} = cell(m_tot, 1);
-            end
+            B = [];
+            if optFlag
+                B = cell(dim, 1);
+                for i = 1:dim
+                    B{i} = cell(m_tot, 1);
+                end
+
+                B0 = sparse(m_tot, m_tot);
+                for i = 1:dim
+                    for j = 1:m_tot
+                        B{i}{j} = B0;
+                    end
+                end
+
+                ind = grid.funcToMatrix(g, 1:m_tot);
 
-            for i = 1:dim
-                for j = 1:m_tot
-                    B{i}{j} = sparse(m_tot, m_tot);
+                % Direction 1
+                for k = 1:m(1)
+                    c = sparse(m(1),1);
+                    c(k) = 1;
+                    [~, B_1D] = ops{1}.D2(c);
+                    for l = 1:m(2)
+                        p = ind(:,l);
+                        B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
+                    end
+                end
+
+                % Direction 2
+                for k = 1:m(2)
+                    c = sparse(m(2),1);
+                    c(k) = 1;
+                    [~, B_1D] = ops{2}.D2(c);
+                    for l = 1:m(1)
+                        p = ind(l,:);
+                        B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
+                    end
                 end
             end
-
-            ind = grid.funcToMatrix(g, 1:m_tot);
+            obj.B = B;
 
-            % Direction 1
-            for k = 1:m(1)
-                c = sparse(m(1),1);
-                c(k) = 1;
-                [~, B_1D] = ops{1}.D2(c);
-                for l = 1:m(2)
-                    p = ind(:,l);
-                    B{1}{(k-1)*m(2) + l}(p, p) = B_1D;
-                end
-            end
-
-            % Direction 2
-            for k = 1:m(2)
-                c = sparse(m(2),1);
-                c(k) = 1;
-                [~, B_1D] = ops{2}.D2(c);
-                for l = 1:m(1)
-                    p = ind(l,:);
-                    B{2}{(l-1)*m(2) + k}(p, p) = B_1D;
-                end
-            end
-
-            obj.B = B;
 
         end