--- a/+multiblock/DiffOp.m	Thu Feb 17 18:55:11 2022 +0100
+++ b/+multiblock/DiffOp.m	Thu Mar 10 16:54:26 2022 +0100
@@ -129,19 +129,20 @@
 
         % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
         function op = getBoundaryOperator(obj, opName, boundary)
+            blockmatrixDiv = obj.blockmatrixDiv{1};
 
             switch class(boundary)
                 case 'cell'
                     blockId = boundary{1};
                     localOp = obj.diffOps{blockId}.getBoundaryOperator(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 = sparse(sum(blockmatrixDiv),0);
                     for i = 1:length(boundary)
                         op = [op, obj.getBoundaryOperator(opName, boundary{i})];
                     end
@@ -208,6 +209,57 @@
             error('not implemented')
         end
 
+        function penalties = interfaceForcing(obj, boundary, neighbour_boundary, type)
+            default_arg('type', []);
+            switch class(boundary)
+                case 'cell'
+                    penalties = obj.singleInterfaceForcing(boundary, neighbour_boundary, type);
+                case 'multiblock.BoundaryGroup'
+                    [n,m] = size(obj.D);
+                    penaltyZero = sparse(n,0);
+
+                    for i = 1:length(boundary)
+                        penaltyParts = obj.interfaceForcing(boundary{i}, neighbour_boundary{i}, type);
+
+                        % Initalize if this is the first boundary
+                        if i==1
+                            penalties = cell(numel(penaltyParts), 1);
+                            for j = 1:numel(penaltyParts)
+                                penalties{j} = penaltyZero;
+                            end
+                        end
+
+                        for j = 1:numel(penaltyParts)
+                            penalties{j} = [penalties{j}, penaltyParts{j}];
+                        end
+
+                    end
+                otherwise
+                    error('Unknown boundary indentifier')
+            end
+
+        end
+
+        function penalties = singleInterfaceForcing(obj, boundary, neighbour_boundary, type)
+            I = boundary{1};
+            b1 = boundary{2};
+
+            J = neighbour_boundary{1};
+            b2 = neighbour_boundary{2};
+
+            % Get local penalties
+            [~, ~, blockPenalties] = obj.diffOps{I}.interface(b1, obj.diffOps{J}, b2, type);
+            % [~, ~, blockPenalties2] = obj.diffOps{J}.interface(b2, b1, obj.diffOps{I});
+
+            % Expand to matrices for full domain.
+            n = numel(blockPenalties);
+            penalties = cell(n, 1);
+            for i = 1:n
+                penalties{i} = multiblock.local2globalPenalty(blockPenalties{i}, obj.blockmatrixDiv, I);
+                             % + multiblock.local2globalPenalty(blockPenalties2{i}, obj.blockmatrixDiv, J);
+            end
+        end
+
         % Size returns the number of degrees of freedom
         function N = size(obj)
             N = 0;