--- a/+multiblock/DiffOp.m	Sun Nov 04 12:36:30 2018 -0800
+++ b/+multiblock/DiffOp.m	Wed Nov 21 18:29:29 2018 -0800
@@ -53,7 +53,11 @@
 
 
             % Build the differentiation matrix
-            obj.blockmatrixDiv = {g.Ns, g.Ns};
+            Ns = zeros(nBlocks,1);
+            for i = 1:nBlocks
+                Ns(i) = length(obj.diffOps{i}.D);
+            end
+            obj.blockmatrixDiv = {Ns, Ns};
             D = blockmatrix.zero(obj.blockmatrixDiv);
             for i = 1:nBlocks
                 D{i,i} = obj.diffOps{i}.D;
@@ -117,7 +121,7 @@
 
         function ops = splitOp(obj, op)
             % Splits a matrix operator into a cell-matrix of matrix operators for
-            % each g.
+            % each grid.
             ops = sparse2cell(op, obj.NNN);
         end
 
@@ -144,6 +148,49 @@
             end
         end
 
+        % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
+        function op = getBoundaryOperatorWrapper(obj, opName, boundary)
+            switch class(boundary)
+                case 'cell'
+                    blockId = boundary{1};
+                    localOp = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2});
+
+                    div = {obj.blockmatrixDiv{1}, size(localOp,2)};
+                    blockOp = blockmatrix.zero(div);
+                    blockOp{blockId,1} = localOp;
+                    op = blockmatrix.toMatrix(blockOp);
+                    return
+                case 'multiblock.BoundaryGroup'
+                    op = sparse(size(obj.D,1),0);
+                    for i = 1:length(boundary)
+                        op = [op, obj.getBoundaryOperatorWrapper(opName, boundary{i})];
+                    end
+                otherwise
+                    error('Unknown boundary indentifier')
+            end
+        end
+
+        function op = getBoundaryQuadrature(obj, boundary)
+            opName = 'H';
+            switch class(boundary)
+                case 'cell'
+                    localOpName = [opName '_' boundary{2}];
+                    blockId = boundary{1};
+                    op = obj.diffOps{blockId}.(localOpName);
+
+                    return
+                case 'multiblock.BoundaryGroup'
+                    N = length(boundary);
+                    H_bm = cell(N,N);
+                    for i = 1:N
+                        H_bm{i,i} = obj.getBoundaryQuadrature(boundary{i});
+                    end
+                    op = blockmatrix.toMatrix(H_bm);
+                otherwise
+                    error('Unknown boundary indentifier')
+            end
+        end
+
         % Creates the closure and penalty matrix for a given boundary condition,
         %    boundary -- the name of the boundary on the form {id,name} where
         %                id is the number of a block and name is the name of a