--- a/+scheme/Elastic2dVariable.m	Wed Oct 24 16:16:43 2018 -0700
+++ b/+scheme/Elastic2dVariable.m	Wed Oct 24 16:17:32 2018 -0700
@@ -43,22 +43,33 @@
         % Kroneckered norms and coefficients
         RHOi_kron
         Hi_kron
+
+        % Structures used for adjoint optimization
+        B
     end
 
     methods
 
-        function obj = Elastic2dVariable(g ,order, lambda_fun, mu_fun, rho_fun, opSet)
+        % The coefficients can either be function handles or grid functions
+        function obj = Elastic2dVariable(g ,order, lambda, mu, rho, opSet)
             default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable});
-            default_arg('lambda_fun', @(x,y) 0*x+1);
-            default_arg('mu_fun', @(x,y) 0*x+1);
-            default_arg('rho_fun', @(x,y) 0*x+1);
+            default_arg('lambda', @(x,y) 0*x+1);
+            default_arg('mu', @(x,y) 0*x+1);
+            default_arg('rho', @(x,y) 0*x+1);
             dim = 2;
 
             assert(isa(g, 'grid.Cartesian'))
 
-            lambda = grid.evalOn(g, lambda_fun);
-            mu = grid.evalOn(g, mu_fun);
-            rho = grid.evalOn(g, rho_fun);
+            if isa(lambda, 'function_handle')
+                lambda = grid.evalOn(g, lambda);
+            end
+            if isa(mu, 'function_handle')
+                mu = grid.evalOn(g, mu);
+            end
+            if isa(rho, 'function_handle')
+                rho = grid.evalOn(g, rho);
+            end
+
             m = g.size();
             m_tot = g.N();
 
@@ -202,7 +213,7 @@
                 end
             end
             obj.D = D;
-            %=========================================%
+            %=========================================%'
 
             % Numerical traction operators for BC.
             % Because d1 =/= e0^T*D1, the numerical tractions are different
@@ -264,6 +275,17 @@
             obj.grid = g;
             obj.dim = dim;
 
+            % Used for adjoint optimization
+            obj.B = cell(1,dim);
+            for i = 1:dim
+                obj.B{i} = zeros(m(i),m(i),m(i));
+                for k = 1:m(i)
+                    c = sparse(m(i),1);
+                    c(k) = 1;
+                    [~, obj.B{i}(:,:,k)] = ops{i}.D2(c);
+                end
+            end
+
         end
 
 
@@ -496,6 +518,28 @@
                             case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
                                 varargout{i} = obj.tau_r{j};
                         end
+                    case 'alpha'
+                        % alpha = alpha(i,j) is the penalty strength for displacement BC. 
+                        tuning = 1.2;
+                        LAMBDA = obj.LAMBDA;
+                        MU = obj.MU;
+
+                        phi = obj.phi{j};
+                        h = obj.h(j);
+                        h11 = obj.H11{j}*h;
+                        gamma = obj.gamma{j};
+                        dim = obj.dim;
+
+                        a_lambda = dim/h11 + 1/(h11*phi);
+                        a_mu_i = 2/(gamma*h);
+                        a_mu_ij = 2/h11 + 1/(h11*phi);
+
+                        d = @kroneckerDelta;  % Kronecker delta
+                        db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
+                        alpha = @(i,k) d(i,k)*tuning*( d(i,j)* a_lambda*LAMBDA ...
+                                                     + d(i,j)* a_mu_i*MU ...
+                                                     + db(i,j)*a_mu_ij*MU );
+                        varargout{i} = alpha;
                     otherwise
                         error(['No such operator: operator = ' op{i}]);
                 end