--- a/+scheme/bcSetup.m	Thu May 24 16:40:12 2018 +0200
+++ b/+scheme/bcSetup.m	Mon Jun 04 17:47:32 2018 +0200
@@ -3,43 +3,70 @@
 % Each bc is a struct with the fields
 %  * type     -- Type of boundary condition
 %  * boundary -- Boundary identifier
-%  * data     -- A function_handle with time and space coordinates as a parameters, for example f(t,x,y) for a 2D problem
+%  * data     -- A function_handle for a function which provides boundary data.(see below)
 % Also takes S_sign which modifies the sign of S, [-1,1]
-% Returns a closure matrix and a forcing function S
+% Returns a closure matrix and a forcing function S.
+%
+% The boundary data function can either be a function of time or a function of time and space coordinates.
+% In the case where it only depends on time it should return the data as grid function for the boundary.
+% In the case where it also takes space coordinates the number of space coordinates should match the number of dimensions of the problem domain.
+% For example in the 2D case: f(t,x,y).
 function [closure, S] = bcSetup(diffOp, bc, S_sign)
     default_arg('S_sign', 1);
     assertType(bc, 'cell');
     assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1');
 
 
+    % Setup storage arrays
     closure = spzeros(size(diffOp));
-    penalties = {};
-    dataFunctions = {};
-    dataParams = {};
+    gridDataPenalties = {};
+    gridDataFunctions = {};
+    symbolicDataPenalties = {};
+    symbolicDataFunctions = {};
+    symbolicDataCoords = {};
 
+    % Collect closures, penalties and data
     for i = 1:length(bc)
         assertType(bc{i}, 'struct');
         [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
         closure = closure + localClosure;
 
         if ~isfield(bc{i},'data') || isempty(bc{i}.data)
+            % Skip to next loop if there is no data
             continue
         end
         assertType(bc{i}.data, 'function_handle');
 
-        coord = diffOp.grid.getBoundary(bc{i}.boundary);
-        assertNumberOfArguments(bc{i}.data, 1+size(coord,2));
+        % Find dimension
+        dim = size(diffOp.grid.getBoundary(bc{i}.boundary), 2);
 
-        penalties{end+1} = penalty;
-        dataFunctions{end+1} = bc{i}.data;
-        dataParams{end+1} = num2cell(coord ,1);
+        if nargin(bc{i}.data) == 1
+            % Grid data
+            boundarySize = [size(diffOp.grid.getBoundary(bc{i}.boundary),1),1];
+            assert_size(bc{i}.data(0), boundarySize); % Eval for t = 0 and make sure the function returns a grid vector of the correct size.
+            gridDataPenalties{end+1} = penalty;
+            gridDataFunctions{end+1} = bc{i}.data;
+        elseif nargin(bc{i}.data) == 1+dim
+            % Symbolic data
+            coord = diffOp.grid.getBoundary(bc{i}.boundary);
+            symbolicDataPenalties{end+1} = penalty;
+            symbolicDataFunctions{end+1} = bc{i}.data;
+            symbolicDataCoords{end+1} = num2cell(coord ,1);
+        else
+            error('sbplib:scheme:bcSetup:DataWrongNumberOfArguments', 'bc{%d}.data has the wrong number of input arguments. Must be either only time or time and space.', i);
+        end
     end
 
+    % Setup penalty function
     O = spzeros(size(diffOp),1);
     function v = S_fun(t)
         v = O;
-        for i = 1:length(dataFunctions)
-            v = v + penalties{i}*dataFunctions{i}(t, dataParams{i}{:});
+        for i = 1:length(gridDataFunctions)
+            v = v + gridDataPenalties{i}*gridDataFunctions{i}(t);
+        end
+
+        for i = 1:length(symbolicDataFunctions)
+            v = v + symbolicDataPenalties{i}*symbolicDataFunctions{i}(t, symbolicDataCoords{i}{:});
         end
 
         v = S_sign * v;