--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/bcSetup.m	Fri Oct 20 16:01:43 2017 +0200
@@ -0,0 +1,48 @@
+% 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)
+    default_arg('S_sign', 1);
+    assertType(bc, '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;
+end