--- a/+scheme/Hypsyst2dCurve.m	Tue Jan 08 11:51:24 2019 +0100
+++ b/+scheme/Hypsyst2dCurve.m	Sat Jan 12 11:57:50 2019 -0800
@@ -169,31 +169,28 @@
             Y = obj.Y;
             xi = obj.xi;
             eta = obj.eta;
+            e_ = obj.getBoundaryOperator('e', boundary);
 
             switch boundary
                 case {'w','W','west'}
-                    e_ = obj.e_w;
                     mat = obj.Ahat;
                     boundPos = 'l';
                     Hi = obj.Hxii;
                     [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_w),Y(obj.index_w),obj.X_eta(obj.index_w),obj.Y_eta(obj.index_w));
                     side = max(length(eta));
                 case {'e','E','east'}
-                    e_ = obj.e_e;
                     mat = obj.Ahat;
                     boundPos = 'r';
                     Hi = obj.Hxii;
                     [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_e),Y(obj.index_e),obj.X_eta(obj.index_e),obj.Y_eta(obj.index_e));
                     side = max(length(eta));
                 case {'s','S','south'}
-                    e_ = obj.e_s;
                     mat = obj.Bhat;
                     boundPos = 'l';
                     Hi = obj.Hetai;
                     [V,Vi,D,signVec] = obj.matrixDiag(mat,X(obj.index_s),Y(obj.index_s),obj.X_xi(obj.index_s),obj.Y_xi(obj.index_s));
                     side = max(length(xi));
                 case {'n','N','north'}
-                    e_ = obj.e_n;
                     mat = obj.Bhat;
                     boundPos = 'r';
                     Hi = obj.Hetai;
@@ -374,5 +371,58 @@
             Vi = [Vi(poseig,:); Vi(zeroeig,:); Vi(negeig,:)];
             signVec = [sum(poseig),sum(zeroeig),sum(negeig)];
         end
+
+        % Returns the boundary operator op for the boundary specified by the string boundary.
+        % op        -- string or a cell array of strings
+        % boundary  -- string
+        function varargout = getBoundaryOperator(obj, op, boundary)
+
+            if ~iscell(op)
+                op = {op};
+            end
+
+            for i = 1:numel(op)
+                switch op{i}
+                case 'e'
+                    switch boundary
+                    case 'w'
+                        e = obj.e_w;
+                    case 'e'
+                        e = obj.e_e;
+                    case 's'
+                        e = obj.e_s;
+                    case 'n'
+                        e = obj.e_n;
+                    otherwise
+                        error('No such boundary: boundary = %s',boundary);
+                    end
+                    varargout{i} = e;
+                end
+            end
+        end
+
+        % Returns square boundary quadrature matrix, of dimension
+        % corresponding to the number of boundary points
+        %
+        % boundary -- string
+        function H_b = getBoundaryQuadrature(obj, boundary)
+
+            e = obj.getBoundaryOperator('e', boundary);
+
+            switch boundary
+                case 'w'
+                    H_b = inv(e'*obj.Hetai*e);
+                case 'e'
+                    H_b = inv(e'*obj.Hetai*e);
+                case 's'
+                    H_b = inv(e'*obj.Hxii*e);
+                case 'n'
+                    H_b = inv(e'*obj.Hxii*e);
+                otherwise
+                    error('No such boundary: boundary = %s',boundary);
+            end
+        end
+
+
     end
 end
\ No newline at end of file