--- a/+scheme/elasticShearVariable.m	Thu Dec 14 13:54:20 2017 -0800
+++ b/+scheme/elasticShearVariable.m	Fri Dec 15 16:23:10 2017 -0800
@@ -4,19 +4,25 @@
         h % Grid spacing
 
         grid
+        dim
 
         order % Order accuracy for the approximation
 
-        a % Variable coefficient lambda of the operator
-        rho % Density
+        A % Variable coefficient lambda of the operator (as diagonal matrix here)
+        RHO % Density (as diagonal matrix here)
 
         D % Total operator
         D1 % First derivatives
         D2 % Second derivatives
         H, Hi % Inner products
         e_l, e_r
-        d_l, d_r % Normal derivatives at the boundary
+        d1_l, d1_r % Normal derivatives at the boundary
+        E % E{i}^T picks out component i
+        
         H_boundary % Boundary inner products
+
+        A_boundary_l % Variable coefficient at boundaries
+        A_boundary_r % 
     end
 
     methods
@@ -56,7 +62,7 @@
             d1_r = cell(dim,1);
 
             for i = 1:dim
-                I{i} = speye{m(i));
+                I{i} = speye(m(i));
                 D1{i} = ops{i}.D1;
                 D2{i} = ops{i}.D2;
                 H{i} =  ops{i}.H;
@@ -69,7 +75,10 @@
 
             %====== Assemble full operators ========
             A = spdiag(a);
+            obj.A = A;
             RHO = spdiag(rho);
+            obj.RHO = RHO;
+
 
             obj.D1 = cell(dim,1);
             obj.D2 = cell(dim,1);
@@ -79,19 +88,19 @@
             obj.d1_r = cell(dim,1);
 
             % D1
-            obj.D1{1} = kron{D1{1},I(2)};
-            obj.D1{2} = kron{I{1},D1(2)};
+            obj.D1{1} = kron(D1{1},I{2});
+            obj.D1{2} = kron(I{1},D1{2});
 
             % Boundary operators
-            obj.e_l{1} = kron{e_l{1},I(2)};
-            obj.e_l{2} = kron{I{1},e_l(2)};
-            obj.e_r{1} = kron{e_r{1},I(2)};
-            obj.e_r{2} = kron{I{1},e_r(2)};
+            obj.e_l{1} = kron(e_l{1},I{2});
+            obj.e_l{2} = kron(I{1},e_l{2});
+            obj.e_r{1} = kron(e_r{1},I{2});
+            obj.e_r{2} = kron(I{1},e_r{2});
 
-            obj.d1_l{1} = kron{d1_l{1},I(2)};
-            obj.d1_l{2} = kron{I{1},d1_l(2)};
-            obj.d1_r{1} = kron{d1_r{1},I(2)};
-            obj.d1_r{2} = kron{I{1},d1_r(2)};
+            obj.d1_l{1} = kron(d1_l{1},I{2});
+            obj.d1_l{2} = kron(I{1},d1_l{2});
+            obj.d1_r{1} = kron(d1_r{1},I{2});
+            obj.d1_r{2} = kron(I{1},d1_r{2});
 
             % D2
             for i = 1:dim
@@ -113,6 +122,7 @@
 
             % Quadratures
             obj.H = kron(H{1},H{2});
+            obj.Hi = inv(obj.H);
             obj.H_boundary = cell(dim,1);
             obj.H_boundary{1} = H{2};
             obj.H_boundary{2} = H{1};
@@ -121,27 +131,29 @@
             obj.A_boundary_l = cell(dim,1);
             obj.A_boundary_r = cell(dim,1);
             for i = 1:dim
-                obj.A_boundary_l{i} = e_l{i}'*A*e_l{i};
-                obj.A_boundary_r{i} = e_r{i}'*A*e_r{i};
+                obj.A_boundary_l{i} = obj.e_l{i}'*A*obj.e_l{i};
+                obj.A_boundary_r{i} = obj.e_r{i}'*A*obj.e_r{i};
             end
 
             % E{i}^T picks out component i.
             E = cell(dim,1);
-            I = speye{mtot,mtot};
+            I = speye(m_tot,m_tot);
             for i = 1:dim
                 e = sparse(dim,1);
                 e(i) = 1;
-                E{i} = kron(e,I)
+                E{i} = kron(I,e);
             end
             obj.E = E;
 
             % Differentiation matrix D (without SAT)
-            D = 0;
+            D2 = obj.D2;
+            D1 = obj.D1;
+            D = sparse(dim*m_tot,dim*m_tot);
             d = @kroneckerDelta;    % Kronecker delta
-            db = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta
+            db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta
             for i = 1:dim
                 for j = 1:dim
-                    D = D + E{i}*inv(rho)*( d(i,j)*D2{i}*E{j}' +...
+                    D = D + E{i}*inv(RHO)*( d(i,j)*D2{i}*E{j}' +...
                                             db(i,j)*D1{j}*A*D1{i}*E{j}' + ...
                                             D2{j}*E{i}' ...
                                           );
@@ -150,16 +162,12 @@
             obj.D = D;
             %=========================================%
 
-            
-
             % Misc.
             obj.m = m;
             obj.h = h;
             obj.order = order;
             obj.grid = g;
-
-            obj.a = a;
-            obj.b = b;
+            obj.dim = dim;
 
             % obj.gamm_u = h_u*ops_u.borrowing.M.d1;
             % obj.gamm_v = h_v*ops_v.borrowing.M.d1;
@@ -178,7 +186,7 @@
             default_arg('parameter', []);
 
             delta = @kroneckerDelta;    % Kronecker delta
-            delta_b = @(i,j) 1-dij(i,j); % Logical not of Kronecker delta
+            delta_b = @(i,j) 1-delta(i,j); % Logical not of Kronecker delta
 
             % j is the coordinate direction of the boundary
             % nj: outward unit normal component. 
@@ -188,16 +196,18 @@
             switch nj
             case 1
                 e = obj.e_r;
-                d = obj.d_r;
+                d = obj.d1_r;
             case -1
                 e = obj.e_l;
-                d = obj.d_l;
+                d = obj.d1_l;
             end
 
             E = obj.E;
             Hi = obj.Hi;
             H_gamma = obj.H_boundary{j};
             A = obj.A;
+            RHO = obj.RHO;
+            D1 = obj.D1;
 
             switch type
                 % Dirichlet boundary condition
@@ -220,18 +230,19 @@
 
                 % Free boundary condition
                 case {'F','f','Free','free'}
-                    closure = 0;
-                    penalty = 0;
+                    closure = sparse(obj.dim*obj.grid.N,obj.dim*obj.grid.N);
+                    penalty = sparse(obj.dim*obj.grid.N,obj.dim*obj.grid.N);
                     % Loop over components
-                    for i = 1:3
-                        closure = closure + E{i}*(-sign)*Hi*e{j}*H_gamma*(...
+                    for i = 1:obj.dim
+                        closure = closure + E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*(...
                                 e{j}'*A*e{j}*d{j}'*E{i}' + ...
-                                delta(i,j)*e{j}'*A*e{i}*d{i}*E{j}' + ...
-                                delta_b(i,j)*A*D1{i}*E{j}' ...
+                                delta(i,j)*e{j}'*A*e{i}*d{i}'*E{j}' + ...
+                                delta_b(i,j)*e{j}'*A*D1{i}*E{j}' ...
                                 );
-                        penalty = penalty - E{i}*(-sign)*Hi*e{j}*H_gamma;
+                        penalty = penalty - E{i}*inv(RHO)*(-nj)*Hi*e{j}*H_gamma*e{j}'*E{j}';
                     end
 
+
                 % Unknown boundary condition
                 otherwise
                     error('No such boundary condition: type = %s',type);
@@ -246,7 +257,7 @@
             error('Interface not implemented');
         end
 
-        % Ruturns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
         function [j, nj] = get_boundary_number(obj, boundary)
 
             switch boundary
@@ -266,6 +277,39 @@
             end
         end
 
+        % Returns the coordinate number and outward normal component for the boundary specified by the string boundary.
+        function [return_op] = get_boundary_operator(obj, op, boundary)
+
+            switch boundary
+                case {'w','W','west','West', 'e', 'E', 'east', 'East'}
+                    j = 1;
+                case {'s','S','south','South', 'n', 'N', 'north', 'North'}
+                    j = 2;
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+
+            switch op
+                case 'e'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.e_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.e_r{j};
+                    end
+                case 'd'
+                    switch boundary
+                        case {'w','W','west','West','s','S','south','South'}
+                            return_op = obj.d_l{j};
+                        case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'}
+                            return_op = obj.d_r{j};
+                    end
+                otherwise
+                    error(['No such operator: operatr = ' op]);
+            end
+
+        end
+
         function N = size(obj)
             N = prod(obj.m);
         end