--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/evalOn.m	Thu Nov 15 16:36:21 2018 -0800
@@ -0,0 +1,40 @@
+% Takes a function and evaluates it on a grid to return a grid function in the
+% form of a (n*k)x1 vector, where n is the number of grid points and k is the
+% number of components of the function.
+%      g -- Grid to evaluate on.
+%   func -- Function to evaluate. May be a function handle or a constant. If
+%           it is a vector value it has to be provided as a column vector,
+function gf = evalOn(g, func)
+    if ~isa(func, 'function_handle')
+        % We should have a constant.
+        assert(size(func,2) == 1,'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector');
+
+        gf = repmat(func,[g.N, 1]);
+        return
+    end
+    % func should now be a function_handle
+    assert(g.D == nargin(func) || nargin(func) < 0,'grid:evalOn:WrongNumberOfInputs', 'The number of inputs of the function must match the dimension of the domain.')
+
+    x = num2cell(g.points(),1);
+    k = numberOfComponents(func, g.D);
+
+    gf = func(x{:});
+    gf = reorderComponents(gf, k);
+end
+
+% Find the number of vector components of func
+function k = numberOfComponents(func, dim)
+    x0 = num2cell(ones(1, dim));
+    f0 = func(x0{:});
+    assert(size(f0,2) == 1, 'grid:evalOn:VectorValuedWrongDim', 'A vector valued function must be given as a column vector');
+    k = length(f0);
+end
+
+% Reorder the components of the function to sit together
+function gf = reorderComponents(a, k)
+    N = length(a)/k;
+    gf = zeros(N*k, 1);
+    for i = 1:k
+        gf(i:k:end) = a((i-1)*N + 1 : i*N);
+    end
+end