--- a/+multiblock/DiffOp.m	Fri Nov 03 10:43:27 2017 +0100
+++ b/+multiblock/DiffOp.m	Fri Nov 03 10:53:15 2017 +0100
@@ -5,12 +5,12 @@
         diffOps
         D
         H
-        
+
         blockmatrixDiv
     end
-    
+
     methods
-        function obj = DiffOp(doHand, grid, order, doParam,timeDep)
+        function obj = DiffOp(doHand, grid, order, doParam)
             %  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, ...)
@@ -25,13 +25,13 @@
             %            doHand(..., doParam{i}{:}) Otherwise doParam is sent as
             %            extra parameters to all doHand: doHand(..., doParam{:})
             default_arg('doParam', [])
-            
+
             [getHand, getParam] = parseInput(doHand, grid, doParam);
-            
+
             nBlocks = grid.nBlocks();
-            
+
             obj.order = order;
-            
+
             % Create the diffOps for each block
             obj.diffOps = cell(1, nBlocks);
             for i = 1:nBlocks
@@ -42,73 +42,48 @@
                 end
                 obj.diffOps{i} = h(grid.grids{i}, order, p{:});
             end
-            
-            
+
+
             % Build the norm matrix
             H = cell(nBlocks, nBlocks);
             for i = 1:nBlocks
                 H{i,i} = obj.diffOps{i}.H;
             end
             obj.H = blockmatrix.toMatrix(H);
-            
-            
+
+
             % Build the differentiation matrix
-            switch timeDep
-                case {'n','no','N','No'}
-                    obj.blockmatrixDiv = {grid.Ns, grid.Ns};
-                    D = blockmatrix.zero(obj.blockmatrixDiv);
-                    for i = 1:nBlocks
-                        D{i,i} = obj.diffOps{i}.D;
+            obj.blockmatrixDiv = {grid.Ns, grid.Ns};
+            D = blockmatrix.zero(obj.blockmatrixDiv);
+            for i = 1:nBlocks
+                D{i,i} = obj.diffOps{i}.D;
+            end
+
+            for i = 1:nBlocks
+                for j = 1:nBlocks
+                    intf = grid.connections{i,j};
+                    if isempty(intf)
+                        continue
                     end
-                    
-                    for i = 1:nBlocks
-                        for j = 1:nBlocks
-                            intf = grid.connections{i,j};
-                            if isempty(intf)
-                                continue
-                            end
-                            
-                            
-                            [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
-                            D{i,i} = D{i,i} + ii;
-                            D{i,j} = D{i,j} + ij;
-                            
-                            [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
-                            D{j,j} = D{j,j} + jj;
-                            D{j,i} = D{j,i} + ji;
-                        end
-                    end
-                    obj.D = blockmatrix.toMatrix(D);
-                case {'y','yes','Y','Yes'}
-                    for i = 1:nBlocks
-                        D{i,i} = @(t)obj.diffOps{i}.D(t);
-                    end
-                    
-                    for i = 1:nBlocks
-                        for j = 1:nBlocks
-                            intf = grid.connections{i,j};
-                            if isempty(intf)
-                                continue
-                            end
-                            
-                            [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
-                            D{i,i} = @(t)D{i,i}(t) + ii(t);
-                            D{i,j} = ij;
-                            
-                            [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
-                            D{j,j} = @(t)D{j,j}(t) + jj(t);
-                            D{j,i} = ji;
-                        end
-                    end
-                obj.D = D;   
+
+
+                    [ii, ij] = obj.diffOps{i}.interface(intf{1}, obj.diffOps{j}, intf{2});
+                    D{i,i} = D{i,i} + ii;
+                    D{i,j} = D{i,j} + ij;
+
+                    [jj, ji] = obj.diffOps{j}.interface(intf{2}, obj.diffOps{i}, intf{1});
+                    D{j,j} = D{j,j} + jj;
+                    D{j,i} = D{j,i} + ji;
+                end
             end
-            
-            
+            obj.D = blockmatrix.toMatrix(D);
+
+
             function [getHand, getParam] = parseInput(doHand, grid, doParam)
                 if ~isa(grid, 'multiblock.Grid')
                     error('multiblock:DiffOp:DiffOp:InvalidGrid', 'Requires a multiblock grid.');
                 end
-                
+
                 if iscell(doHand) && length(doHand) == grid.nBlocks()
                     getHand = @(i)doHand{i};
                 elseif isa(doHand, 'function_handle')
@@ -116,62 +91,89 @@
                 else
                     error('multiblock:DiffOp:DiffOp:InvalidGridDoHand', 'doHand must be a function handle or a cell array of length grid.nBlocks');
                 end
-                
+
                 if isempty(doParam)
                     getParam = @(i){};
                     return
                 end
-                
+
                 if ~iscell(doParam)
                     getParam = @(i)doParam;
                     return
                 end
-                
+
                 % doParam is a non-empty cell-array
-                
+
                 if length(doParam) == grid.nBlocks() && all(cellfun(@iscell, doParam))
                     % doParam is a cell-array of cell-arrays
                     getParam = @(i)doParam{i};
                     return
                 end
-                
+
                 getParam = @(i)doParam;
             end
         end
-        
+
         function ops = splitOp(obj, op)
             % Splits a matrix operator into a cell-matrix of matrix operators for
             % each grid.
             ops = sparse2cell(op, obj.NNN);
         end
-        
-        function op = getBoundaryOperator(obj, op, boundary)
-            if iscell(boundary)
-                localOpName = [op '_' boundary{2}];
-                blockId = boundary{1};
-                localOp = obj.diffOps{blockId}.(localOpName);
-                
-                div = {obj.blockmatrixDiv{1}, size(localOp,2)};
-                blockOp = blockmatrix.zero(div);
-                blockOp{blockId,1} = localOp;
-                op = blockmatrix.toMatrix(blockOp);
-                return
-            else
-                % Boundary är en sträng med en boundary group i.
+
+        % Get a boundary operator specified by opName for the given boundary/BoundaryGroup
+        function op = getBoundaryOperator(obj, opName, boundary)
+            switch class(boundary)
+                case 'cell'
+                    localOpName = [opName '_' boundary{2}];
+                    blockId = boundary{1};
+                    localOp = obj.diffOps{blockId}.(localOpName);
+
+                    div = {obj.blockmatrixDiv{1}, 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);
+                    for i = 1:length(boundary)
+                        op = [op, obj.getBoundaryOperator(opName, boundary{i})];
+                    end
+                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
-        %                boundary of that block example: {1,'s'} or {3,'w'}
+        %                boundary of that block example: {1,'s'} or {3,'w'}. It
+        %                can also be a boundary group
         function [closure, penalty] = boundary_condition(obj, boundary, type)
+            switch class(boundary)
+                case 'cell'
+                    [closure, penalty] = obj.singleBoundaryCondition(boundary, type);
+                case 'multiblock.BoundaryGroup'
+                    [n,m] = size(obj.D);
+                    closure = sparse(n,m);
+                    penalty = sparse(n,0);
+                    for i = 1:length(boundary)
+                        [closurePart, penaltyPart] = obj.boundary_condition(boundary{i}, type);
+                        closure = closure + closurePart;
+                        penalty = [penalty, penaltyPart];
+                    end
+                otherwise
+                    error('Unknown boundary indentifier')
+            end
+
+        end
+
+        function [closure, penalty] = singleBoundaryCondition(obj, boundary, type)
             I = boundary{1};
             name = boundary{2};
-            
+
             % Get the closure and penaly matrices
             [blockClosure, blockPenalty] = obj.diffOps{I}.boundary_condition(name, type);
-            
+
             % Expand to matrix for full domain.
             div = obj.blockmatrixDiv;
             if ~iscell(blockClosure)
@@ -185,25 +187,27 @@
                     closure{i} = blockmatrix.toMatrix(temp);
                 end
             end
-            
-            div{2} = size(blockPenalty, 2); % Penalty is a column vector
+
             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
         end
-        
+
         function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary)
-            
+            error('not implemented')
         end
-        
+
         % Size returns the number of degrees of freedom
         function N = size(obj)
             N = 0;