--- a/+scheme/bcSetup.m	Tue Dec 04 14:54:28 2018 -0800
+++ b/+scheme/bcSetup.m	Tue Dec 04 15:24:36 2018 -0800
@@ -1,48 +1,20 @@
-% function [closure, S] = bcSetup(diffOp, bc)
 % Takes a diffOp and a cell array of boundary condition definitions.
 % 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
-% Also takes S_sign which modifies the sign of S, [-1,1]
-% Returns a closure matrix and a forcing function S
-function [closure, S] = bcSetup(diffOp, bc, S_sign)
+%  * data     -- A function_handle for a function which provides boundary data.(see below)
+% Also takes S_sign which modifies the sign of the penalty function, [-1,1]
+% 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, bcs, S_sign)
     default_arg('S_sign', 1);
-    assertType(bc, 'cell');
+    assertType(bcs, 'cell');
     assert(S_sign == 1 || S_sign == -1, 'S_sign must be either 1 or -1');
 
-
-    closure = spzeros(size(diffOp));
-    penalties = {};
-    dataFunctions = {};
-    dataParams = {};
-
-    for i = 1:length(bc)
-        assertType(bc{i}, 'struct');
-        [localClosure, penalty] = diffOp.boundary_condition(bc{i}.boundary, bc{i}.type);
-        closure = closure + localClosure;
-
-        if isempty(bc{i}.data)
-            continue
-        end
-        assertType(bc{i}.data, 'function_handle');
-
-        coord = diffOp.grid.getBoundary(bc{i}.boundary);
-        assertNumberOfArguments(bc{i}.data, 1+size(coord,2));
-
-        penalties{end+1} = penalty;
-        dataFunctions{end+1} = bc{i}.data;
-        dataParams{end+1} = num2cell(coord ,1);
-    end
-
-    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}{:});
-        end
-
-        v = S_sign * v;
-    end
-    S = @S_fun;
+    [closure, penalties] = scheme.bc.closureSetup(diffOp, bcs);
+    S = scheme.bc.forcingSetup(diffOp, penalties, bcs, S_sign);
 end