diff -r 2b03ee11e48b -r 31abb4b11133 +grid/TODO.txt
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/TODO.txt	Fri Oct 20 23:31:20 2017 +0200
@@ -0,0 +1,1 @@
+% TODO: Rename grid package. name conflicts with built in function

diff -r 2b03ee11e48b -r 31abb4b11133 +multiblock/DiffOp.m
--- a/+multiblock/DiffOp.m	Fri Oct 20 23:30:05 2017 +0200
+++ b/+multiblock/DiffOp.m	Fri Oct 20 23:31:20 2017 +0200
@@ -194,6 +194,7 @@
                 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);

diff -r 2b03ee11e48b -r 31abb4b11133 +noname/Animation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+noname/Animation.m	Fri Oct 20 23:31:20 2017 +0200
@@ -0,0 +1,75 @@
+classdef Animation < handle
+    properties
+        timeStepper
+        representationMaker
+        updaters
+    end
+
+    % add input validation
+
+    methods
+        function obj = Animation(timeStepper, representationMaker, updaters);
+            obj.timeStepper = timeStepper;
+            obj.updaters = updaters;
+            obj.representationMaker = representationMaker;
+        end
+
+        function update(obj, r)
+            for i = 1:length(obj.updaters)
+                obj.updaters{i}(r);
+            end
+            drawnow
+        end
+
+        function run(obj, tEnd, timeModifier, do_pause)
+            default_arg('do_pause', false)
+
+            function next_t = G(next_t)
+                obj.timeStepper.evolve(next_t);
+                r = obj.representationMaker(obj.timeStepper);
+                obj.update(r);
+
+                if do_pause
+                    pause
+                end
+            end
+
+            anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier);
+        end
+
+        function step(obj, tEnd, do_pause)
+            default_arg('do_pause', false)
+
+            while obj.timeStepper.t < tEnd
+                obj.timeStepper.step();
+
+                r = obj.representationMaker(obj.timeStepper);
+                obj.update(r);
+
+                % TODO: Make it never go faster than a certain fram rate
+
+                if do_pause
+                    pause
+                end
+            end
+        end
+
+        function saveMovie(obj, tEnd, timeModifier, figureHandle, dirname)
+            save_frame = anim.setup_fig_mov(figureHandle, dirname);
+
+            function next_t = G(next_t)
+                obj.timeStepper.evolve(next_t);
+                r = obj.representationMaker(obj.timeStepper);
+                obj.update(r);
+
+                save_frame();
+            end
+
+            fprintf('Generating and saving frames to: ..\n')
+            anim.animate(@G, obj.timeStepper.t, tEnd, timeModifier);
+            fprintf('Generating movies...\n')
+            cmd = sprintf('bash %s/+anim/make_movie.sh %s', sbplibLocation(),dirname);
+            system(cmd);
+        end
+    end
+end

diff -r 2b03ee11e48b -r 31abb4b11133 +noname/calculateErrors.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+noname/calculateErrors.m	Fri Oct 20 23:31:20 2017 +0200
@@ -0,0 +1,40 @@
+% [discr, trueSolution] =  schemeFactory(m)
+%     where trueSolution should be a timeSnapshot of the true solution a time T
+% T is the end time
+% m are grid size parameters.
+% N are number of timesteps to use for each gird size
+% timeOpt are options for the timeStepper
+function e = calculateErrors(schemeFactory, T, m, N, errorFun, timeOpt)
+    assertType(schemeFactory, 'function_handle');
+    assertNumberOfArguments(schemeFactory, 1);
+    assertScalar(T);
+    assert(length(m) == length(N), 'Vectors m and N must have the same length');
+    assertType(errorFun, 'function_handle');
+    assertNumberOfArguments(errorFun, 2);
+    default_arg('timeOpt');
+
+    e = [];
+    for i = 1:length(m)
+        done = timeTask('m = %3d ', m(i));
+
+        [discr, trueSolution] = schemeFactory(m(i));
+
+        timeOpt.k = T/N(i);
+        ts = discr.getTimestepper(timeOpt);
+        ts.stepTo(N(i), true);
+        approxSolution = discr.getTimeSnapshot(ts);
+
+        e(i) = errorFun(trueSolution, approxSolution);
+
+        fprintf('e = %.4e', e(i))
+        done()
+    end
+    fprintf('\n')
+end
+
+
+%% Example error function
+% u_true = grid.evalOn(dr.grid, @(x,y)trueSolution(T,x,y));
+% err = u_true-u_false;
+% e(i) = norm(err)/norm(u_true);
+% % e(i) = sqrt(err'*d.H*d.J*err/(u_true'*d.H*d.J*u_true));

diff -r 2b03ee11e48b -r 31abb4b11133 +scheme/Beam.m
--- a/+scheme/Beam.m	Fri Oct 20 23:30:05 2017 +0200
+++ b/+scheme/Beam.m	Fri Oct 20 23:31:20 2017 +0200
@@ -126,6 +126,44 @@
                     penalty{1} = -obj.Hi*tau;
                     penalty{1} = -obj.Hi*sig;
 
+                case 'e'
+                    alpha = obj.alpha;
+                    tuning = 1.1;
+
+                    tau1 = tuning * alpha/delt;
+                    tau4 = s*alpha;
+
+                    tau = tau1*e+tau4*d3;
+
+                    closure = obj.Hi*tau*e';
+                    penalty = -obj.Hi*tau;
+                case 'd1'
+                    alpha = obj.alpha;
+
+                    tuning = 1.1;
+
+                    sig2 = tuning * alpha/gamm;
+                    sig3 = -s*alpha;
+
+                    sig = sig2*d1+sig3*d2;
+
+                    closure = obj.Hi*sig*d1';
+                    penalty = -obj.Hi*sig;
+
+                case 'd2'
+                    a = obj.alpha;
+
+                    tau =  s*a*d1;
+
+                    closure = obj.Hi*tau*d2';
+                    penalty = -obj.Hi*tau;
+                case 'd3'
+                    a = obj.alpha;
+
+                    sig = -s*a*e;
+
+                    closure = obj.Hi*sig*d3';
+                    penalty = -obj.Hi*sig;
 
                 otherwise % Unknown, boundary condition
                     error('No such boundary condition: type = %s',type);

diff -r 2b03ee11e48b -r 31abb4b11133 +scheme/TODO.txt
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/TODO.txt	Fri Oct 20 23:31:20 2017 +0200
@@ -0,0 +1,1 @@
+% TODO: Rename package and abstract class to diffOp

diff -r 2b03ee11e48b -r 31abb4b11133 +scheme/bcSetup.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+scheme/bcSetup.m	Fri Oct 20 23:31:20 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

diff -r 2b03ee11e48b -r 31abb4b11133 +time/Timestepper.m
--- a/+time/Timestepper.m	Fri Oct 20 23:30:05 2017 +0200
+++ b/+time/Timestepper.m	Fri Oct 20 23:31:20 2017 +0200
@@ -62,6 +62,7 @@
 
 
         function [v, t] = stepTo(obj, n, progress_bar)
+            assertScalar(n);
             default_arg('progress_bar',false);
 
             [v, t] = obj.stepN(n-obj.n, progress_bar);

diff -r 2b03ee11e48b -r 31abb4b11133 assertScalar.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/assertScalar.m	Fri Oct 20 23:31:20 2017 +0200
@@ -0,0 +1,5 @@
+function assertScalar(obj)
+    if ~isscalar(obj)
+        error('sbplib:assertScalar:notScalar', '"%s" must be scalar, found size "%s"', inputname(1), toString(size(obj)));
+    end
+end

diff -r 2b03ee11e48b -r 31abb4b11133 diffSymfun.m
--- a/diffSymfun.m	Fri Oct 20 23:30:05 2017 +0200
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,7 +0,0 @@
-% Differentiates a symbolic function like diff does, but keeps the function as a symfun
-function g = diffSymfun(f, varargin)
-    assertType(f, 'symfun');
-
-    args = argnames(f);
-    g = symfun(diff(f,varargin{:}), args);
-end

diff -r 2b03ee11e48b -r 31abb4b11133 sbplibLocation.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/sbplibLocation.m	Fri Oct 20 23:31:20 2017 +0200
@@ -0,0 +1,4 @@
+function location = sbplibLocation()
+    scriptname  = mfilename('fullpath');
+    [location, ~, ~] = fileparts(scriptname);
+end

diff -r 2b03ee11e48b -r 31abb4b11133 sbplibVersion.m
--- a/sbplibVersion.m	Fri Oct 20 23:30:05 2017 +0200
+++ b/sbplibVersion.m	Fri Oct 20 23:31:20 2017 +0200
@@ -1,11 +1,10 @@
 % Prints the version and location of the sbplib currently in use.
 function sbplibVersion()
-    scriptname  = mfilename('fullpath');
-    [folder,~,~] = fileparts(scriptname);
+    location = sbplibLocation();
 
     name = 'sbplib (feature/grids)';
     ver = '0.0.x';
 
     fprintf('%s %s\n', name, ver);
-    fprintf('Running in:\n%s\n',folder);
+    fprintf('Running in:\n%s\n', location);
 end
\ No newline at end of file