--- a/+scheme/LaplaceCurvilinear.m	Tue Aug 29 13:04:07 2017 +0200
+++ b/+scheme/LaplaceCurvilinear.m	Tue Aug 29 13:10:57 2017 +0200
@@ -223,19 +223,13 @@
             default_arg('type','neumann');
             default_arg('parameter', []);
 
-            [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv  ,              ~,          ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
+            [e, d, s, gamm, halfnorm_inv  ,              ~,          ~, ~, scale_factor] = obj.get_boundary_ops(boundary);
             switch type
                 % Dirichlet boundary condition
                 case {'D','d','dirichlet'}
-                    % v denotes the solution in the neighbour domain
                     tuning = 1.2;
                     % tuning = 20.2;
-                    [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
-
-                    A_n = spdiag(a_n);
-                    A_t = spdiag(a_t);
-
-                    F = s*(A_n*d_n' + A_t*d_t')';
+                    [e, F, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t] = obj.get_boundary_ops(boundary);
 
                     u = obj;
 
@@ -254,10 +248,6 @@
 
                 % Neumann boundary condition
                 case {'N','n','neumann'}
-                    A_n = spdiag(a_n);
-                    A_t = spdiag(a_t);
-                    d = s*(A_n * d_n' + A_t*d_t')';
-
                     tau1 = -1;
                     tau2 = 0;
                     tau = s*obj.a*obj.Ji*(tau1*e + tau2*d);
@@ -270,10 +260,6 @@
                     default_arg('parameter', 1);
                     beta = parameter;
 
-                    A_n = spdiag(a_n);
-                    A_t = spdiag(a_t);
-                    d = s*(A_n * d_n' + A_t*d_t')'; % outward facing normal derivative
-
                     tau = -obj.a * 1/beta*obj.Ji*e;
 
                     closure{1} = halfnorm_inv*tau*spdiag(scale_factor)*e';
@@ -291,16 +277,8 @@
             % v denotes the solution in the neighbour domain
             tuning = 1.2;
             % tuning = 20.2;
-            [e_u, d_n_u, d_t_u, a_n_u, a_t_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary);
-            [e_v, d_n_v, d_t_v, a_n_v, a_t_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
-
-            A_n_u = spdiag(a_n_u);
-            A_t_u = spdiag(a_t_u);
-            A_n_v = spdiag(a_n_v);
-            A_t_v = spdiag(a_t_v);
-
-            F_u = s_u*(A_n_u * d_n_u' + A_t_u*d_t_u')';
-            F_v = s_v*(A_n_v * d_n_v' + A_t_v*d_t_v')';
+            [e_u, F_u, s_u, gamm_u, halfnorm_inv_u_n, halfnorm_inv_u_t, halfnorm_u_t, I_u] = obj.get_boundary_ops(boundary);
+            [e_v, F_v, s_v, gamm_v, halfnorm_inv_v_n, halfnorm_inv_v_t, halfnorm_v_t, I_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
 
             u = obj;
             v = neighbour_scheme;
@@ -327,7 +305,7 @@
         % The right boundary is considered the positive boundary
         %
         %  I -- the indecies of the boundary points in the grid matrix
-        function [e, d_n, d_t, a_n, a_t, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
+        function [e, d, s, gamm, halfnorm_inv_n, halfnorm_inv_t, halfnorm_t, I, scale_factor] = get_boundary_ops(obj, boundary)
 
             % gridMatrix = zeros(obj.m(2),obj.m(1));
             % gridMatrix(:) = 1:numel(gridMatrix);
@@ -337,43 +315,31 @@
             switch boundary
                 case 'w'
                     e = obj.e_w;
-                    d_n = obj.du_w;
-                    d_t = obj.dv_w;
+                    d = obj.d_w;
                     s = -1;
 
                     I = ind(1,:);
-                    a_n = obj.a11(I);
-                    a_t = obj.a12(I);
                     scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
                 case 'e'
                     e = obj.e_e;
-                    d_n = obj.du_e;
-                    d_t = obj.dv_e;
+                    d = obj.d_e;
                     s = 1;
 
                     I = ind(end,:);
-                    a_n = obj.a11(I);
-                    a_t = obj.a12(I);
                     scale_factor = sqrt(obj.x_v(I).^2 + obj.y_v(I).^2);
                 case 's'
                     e = obj.e_s;
-                    d_n = obj.dv_s;
-                    d_t = obj.du_s;
+                    d = obj.d_s;
                     s = -1;
 
                     I = ind(:,1)';
-                    a_n = obj.a22(I);
-                    a_t = obj.a12(I);
                     scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
                 case 'n'
                     e = obj.e_n;
-                    d_n = obj.dv_n;
-                    d_t = obj.du_n;
+                    d = obj.d_n;
                     s = 1;
 
                     I = ind(:,end)';
-                    a_n = obj.a22(I);
-                    a_t = obj.a12(I);
                     scale_factor = sqrt(obj.x_u(I).^2 + obj.y_u(I).^2);
                 otherwise
                     error('No such boundary: boundary = %s',boundary);