--- a/+scheme/LaplaceCurvilinear.m	Fri Mar 29 14:50:50 2019 -0700
+++ b/+scheme/LaplaceCurvilinear.m	Wed Apr 10 11:00:27 2019 -0700
@@ -278,12 +278,12 @@
             default_arg('type','neumann');
             default_arg('parameter', []);
 
-            e = obj.getBoundaryOperator('e', boundary);
-            d = obj.getBoundaryOperator('d', boundary);
-            H_b = obj.getBoundaryQuadrature(boundary);
-            s_b = obj.getBoundaryScaling(boundary);
+            e               = obj.getBoundaryOperator('e', boundary);
+            d               = obj.getBoundaryOperator('d', boundary);
+            H_b             = obj.getBoundaryQuadrature(boundary);
+            s_b             = obj.getBoundaryScaling(boundary);
             [th_H, ~, th_R] = obj.getBoundaryBorrowing(boundary);
-            m = obj.getBoundaryNumber(boundary);
+            m               = obj.getBoundaryNumber(boundary);
 
             K = obj.K;
             J = obj.J;
@@ -358,10 +358,10 @@
             v = neighbour_scheme;
 
             % Boundary operators, u
-            e_u = u.getBoundaryOperator('e', boundary);
-            d_u = u.getBoundaryOperator('d', boundary);
-            gamm_u = u.getBoundaryBorrowing(boundary);
-            s_b_u = u.getBoundaryScaling(boundary);
+            e_u     = u.getBoundaryOperator('e', boundary);
+            d_u     = u.getBoundaryOperator('d', boundary);
+            gamm_u  = u.getBoundaryBorrowing(boundary);
+            s_b_u   = u.getBoundaryScaling(boundary);
             [th_H_u, ~, th_R_u] = u.getBoundaryBorrowing(boundary);
             m_u = u.getBoundaryNumber(boundary);
 
@@ -371,10 +371,10 @@
             b_b_u = e_u'*u.b*e_u;
 
             % Boundary operators, v
-            e_v = v.getBoundaryOperator('e', neighbour_boundary);
-            d_v = v.getBoundaryOperator('d', neighbour_boundary);
-            gamm_v = v.getBoundaryBorrowing(neighbour_boundary);
-            s_b_v = v.getBoundaryScaling(neighbour_boundary);
+            e_v     = v.getBoundaryOperator('e', neighbour_boundary);
+            d_v     = v.getBoundaryOperator('d', neighbour_boundary);
+            gamm_v  = v.getBoundaryBorrowing(neighbour_boundary);
+            s_b_v   = v.getBoundaryScaling(neighbour_boundary);
             [th_H_v, ~, th_R_v] = v.getBoundaryBorrowing(neighbour_boundary);
             m_v = v.getBoundaryNumber(neighbour_boundary);
 
@@ -431,16 +431,16 @@
 
             % u denotes the solution in the own domain
             % v denotes the solution in the neighbour domain
-            e_u = obj.getBoundaryOperator('e', boundary);
-            d_u = obj.getBoundaryOperator('d', boundary);
-            H_b_u = obj.getBoundaryQuadrature(boundary);
-            I_u = obj.getBoundaryIndices(boundary);
+            e_u    = obj.getBoundaryOperator('e', boundary);
+            d_u    = obj.getBoundaryOperator('d', boundary);
+            H_b_u  = obj.getBoundaryQuadrature(boundary);
+            I_u    = obj.getBoundaryIndices(boundary);
             gamm_u = obj.getBoundaryBorrowing(boundary);
 
-            e_v = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
-            d_v = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
-            H_b_v = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
-            I_v = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
+            e_v    = neighbour_scheme.getBoundaryOperator('e', neighbour_boundary);
+            d_v    = neighbour_scheme.getBoundaryOperator('d', neighbour_boundary);
+            H_b_v  = neighbour_scheme.getBoundaryQuadrature(neighbour_boundary);
+            I_v    = neighbour_scheme.getBoundaryIndices(neighbour_boundary);
             gamm_v = neighbour_scheme.getBoundaryBorrowing(neighbour_boundary);
 
 
@@ -530,6 +530,8 @@
         % Returns the indices of the boundary points in the grid matrix
         % boundary -- string
         function I = getBoundaryIndices(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
             ind = grid.funcToMatrix(obj.grid, 1:prod(obj.m));
             switch boundary
                 case 'w'
@@ -540,14 +542,14 @@
                     I = ind(:,1)';
                 case 'n'
                     I = ind(:,end)';
-                otherwise
-                    error('No such boundary: boundary = %s',boundary);
             end
         end
 
         % Returns borrowing constant gamma
         % boundary -- string
         function [theta_H, theta_M, theta_R] = getBoundaryBorrowing(obj, boundary)
+            assertIsMember(boundary, {'w', 'e', 's', 'n'})
+
             switch boundary
                 case {'w','e'}
                     theta_H = obj.theta_H_u;
@@ -557,8 +559,6 @@
                     theta_H = obj.theta_H_v;
                     theta_M = obj.theta_M_v;
                     theta_R = obj.theta_R_v;
-                otherwise
-                    error('No such boundary: boundary = %s',boundary);
             end
         end