--- a/+multiblock/DiffOp.m	Mon Dec 17 21:07:59 2018 -0800
+++ b/+multiblock/DiffOp.m	Wed Dec 19 06:54:47 2018 +0100
@@ -149,21 +149,23 @@
         end
 
         % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
-        function op = getBoundaryOperatorWrapper(obj, opName, boundary)
+        function op = getBoundaryOperatorWrapper(obj, opName, boundary, blockmatrixDiv)
+            default_arg('blockmatrixDiv', obj.blockmatrixDiv);
+
             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)};
+                    div = {blockmatrixDiv, 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);
+                    op = [];
                     for i = 1:length(boundary)
-                        op = [op, obj.getBoundaryOperatorWrapper(opName, boundary{i})];
+                        op = [op, obj.getBoundaryOperatorWrapper(opName, boundary{i}, blockmatrixDiv)];
                     end
                 otherwise
                     error('Unknown boundary indentifier')
@@ -195,7 +197,7 @@
                         localOpCell = obj.diffOps{blockId}.get_boundary_operator(opName, boundary{2});
                         localOp = localOpCell{i};
 
-                        div = {blockmatrixDiv, size(localOp,2)}
+                        div = {blockmatrixDiv, size(localOp,2)};
                         blockOp = blockmatrix.zero(div);
                         blockOp{blockId,1} = localOp;
                         op = blockmatrix.toMatrix(blockOp);