--- a/+multiblock/DiffOp.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+multiblock/DiffOp.m	Tue Dec 04 15:24:36 2018 -0800
@@ -10,13 +10,13 @@
     end
 
     methods
-        function obj = DiffOp(doHand, grid, order, doParam, intfTypes)
+        function obj = DiffOp(doHand, g, order, doParam, intfTypes)
             %  doHand -- may either be a function handle or a cell array of
             %            function handles for each grid. The function handle(s)
             %            should be on the form do = doHand(grid, order, ...)
             %            Additional parameters for each doHand may be provided in
             %            the doParam input.
-            %    grid -- a multiblock grid
+            %       g -- a multiblock grid
             %   order -- integer specifying the order of accuracy
             % doParam -- may either be a cell array or a cell array of cell arrays
             %            for each block. If it is a cell array with length equal
@@ -28,12 +28,12 @@
             % intfTypes (optional) -- nBlocks x nBlocks cell array of types for
             %                                 every interface.
             default_arg('doParam', [])
-            default_arg('intfTypes', cell(grid.nBlocks(), grid.nBlocks()) );
+            default_arg('intfTypes', cell(g.nBlocks(), g.nBlocks()) );
 
-            [getHand, getParam] = parseInput(doHand, grid, doParam);
+            [getHand, getParam] = parseInput(doHand, g, doParam);
 
             obj.order = order;
-            nBlocks = grid.nBlocks();
+            nBlocks = g.nBlocks();
 
             % Create the diffOps for each block
             obj.diffOps = cell(1, nBlocks);
@@ -43,7 +43,7 @@
                 if ~iscell(p)
                     p = {p};
                 end
-                obj.diffOps{i} = h(grid.grids{i}, order, p{:});
+                obj.diffOps{i} = h(g.grids{i}, order, p{:});
             end
 
 
@@ -56,7 +56,11 @@
 
 
             % Build the differentiation matrix
-            obj.blockmatrixDiv = {grid.Ns, grid.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;
@@ -64,7 +68,7 @@
 
             for i = 1:nBlocks
                 for j = 1:nBlocks
-                    intf = grid.connections{i,j};
+                    intf = g.connections{i,j};
                     if isempty(intf)
                         continue
                     end
@@ -79,14 +83,15 @@
                 end
             end
             obj.D = blockmatrix.toMatrix(D);
+            obj.grid = g;
 
 
-            function [getHand, getParam] = parseInput(doHand, grid, doParam)
-                if ~isa(grid, 'multiblock.Grid')
+            function [getHand, getParam] = parseInput(doHand, g, doParam)
+                if ~isa(g, 'multiblock.Grid')
                     error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
                 end
 
-                if iscell(doHand) && length(doHand) == grid.nBlocks()
+                if iscell(doHand) && length(doHand) == g.nBlocks()
                     getHand = @(i)doHand{i};
                 elseif isa(doHand, 'function_handle')
                     getHand = @(i)doHand;
@@ -106,7 +111,7 @@
 
                 % doParam is a non-empty cell-array
 
-                if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam))
+                if length(doParam) == g.nBlocks() && all(cellfun(@iscell, doParam))
                     % doParam is a cell-array of cell-arrays
                     getParam = @(i)doParam{i};
                     return
@@ -145,6 +150,27 @@
             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
@@ -177,33 +203,8 @@
             [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
 
             % Expand to matrix for full domain.
-            div = obj.blockmatrixDiv;
-            if ~iscell(blockClosure)
-                temp = blockmatrix.zero(div);
-                temp{I,I} = blockClosure;
-                closure = blockmatrix.toMatrix(temp);
-            else
-                for i = 1:length(blockClosure)
-                    temp = blockmatrix.zero(div);
-                    temp{I,I} = blockClosure{i};
-                    closure{i} = blockmatrix.toMatrix(temp);
-                end
-            end
-
-            if ~iscell(blockPenalty)
-                div{2} = size(blockPenalty, 2); % Penalty is a column vector
-                p = blockmatrix.zero(div);
-                p{I} = blockPenalty;
-                penalty = blockmatrix.toMatrix(p);
-            else
-                % TODO: used by beam equation, should be eliminated. SHould only set one BC per call
-                for i = 1:length(blockPenalty)
-                    div{2} = size(blockPenalty{i}, 2); % Penalty is a column vector
-                    p = blockmatrix.zero(div);
-                    p{I} = blockPenalty{i};
-                    penalty{i} = blockmatrix.toMatrix(p);
-                end
-            end
+            closure = multiblock.local2globalClosure(blockClosure, obj.blockmatrixDiv, I);
+            penalty = multiblock.local2globalPenalty(blockPenalty, obj.blockmatrixDiv, I);
         end
 
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)