--- a/+multiblock/DiffOp.m	Wed Jun 20 12:59:07 2018 +0200
+++ b/+multiblock/DiffOp.m	Wed Jun 20 13:00:01 2018 +0200
@@ -144,6 +144,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

--- a/+multiblock/Laplace.m	Wed Jun 20 12:59:07 2018 +0200
+++ b/+multiblock/Laplace.m	Wed Jun 20 13:00:01 2018 +0200
@@ -42,6 +42,10 @@
             op = getBoundaryOperator(obj.mbDiffOp, opName, boundary);
         end
 
+        function op = getBoundaryQuadrature(obj, boundary)
+            op = getBoundaryQuadrature(obj.mbDiffOp, boundary);
+        end
+
         function [closure, penalty] = boundary_condition(obj,boundary,type) % TODO: Change name to boundaryCondition
             [closure, penalty] = boundary_condition(obj.mbDiffOp, boundary, type);
         end