--- a/+multiblock/DiffOp.m	Sun Jul 26 17:38:28 2020 -0700
+++ b/+multiblock/DiffOp.m	Sun Jul 26 20:05:10 2020 -0700
@@ -209,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;