diff +time/SBPInTimeSecondOrderFormImplicit.m @ 460:e0caae9ef6ed feature/grids

Add SBPinTime for linear DAE formulations (BUGS!)
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 19 Jun 2017 16:50:13 +0200
parents +time/SBPInTimeSecondOrderForm.m@c92d2f8319c2
children 45a3924140f4
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+time/SBPInTimeSecondOrderFormImplicit.m	Mon Jun 19 16:50:13 2017 +0200
@@ -0,0 +1,80 @@
+classdef SBPInTimeSecondOrderFormImplicit < time.Timestepper
+    properties
+        A, B, C, f
+        AA, BB, ff
+
+        n
+        t
+        k
+
+        firstOrderTimeStepper
+    end
+
+    methods
+        % Solves A*u_tt + B*u_t + C*u + f(t) = 0
+        % A, B can either both be constants or both be function handles,
+        % They can also be omitted by setting them equal to the empty matrix.
+        function obj = SBPInTimeSecondOrderFormImplicit(A, B, C, f, k, t0, v0, v0t, TYPE, order, blockSize)
+            default_arg('f', []);
+            default_arg('TYPE', []);
+            default_arg('order', []);
+            default_arg('blockSize',[]);
+
+            m = length(v0);
+
+            default_arg('A', sparse(m, m));
+            default_arg('B', sparse(m, m));
+            default_arg('C', sparse(m, m));
+
+            I = speye(m);
+            O = sparse(m,m);
+
+            % Rewrite to
+            % AA*w_t = BB*w + ff(t);
+
+            obj.AA = [
+                 I, O;
+                 O, A;
+            ];
+            obj.BB = [
+                 O,  I;
+                -B, -C;
+            ];
+
+            if ~isempty(f)
+                obj.ff = @(t)[
+                    sparse(m,1);
+                           -f(t);
+                ];
+            else
+                obj.ff = @(t) sparse(2*m,1);
+            end
+
+            w0 = [v0; v0t];
+
+            obj.k = k;
+            obj.t = t0;
+            obj.n = 0;
+
+            obj.firstOrderTimeStepper = time.SBPInTimeImplicitFormulation(obj.AA, obj.BB, obj.ff, obj.k, obj.t, w0, TYPE, order, blockSize);
+        end
+
+        function [v,t] = getV(obj)
+            w = obj.firstOrderTimeStepper.getV();
+            v = w(1:end/2);
+            t = obj.t;
+        end
+
+        function [vt,t] = getVt(obj)
+            w = obj.firstOrderTimeStepper.getV();
+            vt = w(end/2+1:end);
+            t = obj.t;
+        end
+
+        function obj = step(obj)
+            obj.firstOrderTimeStepper.step();
+            obj.t = obj.firstOrderTimeStepper.t;
+            obj.n = obj.firstOrderTimeStepper.n;
+        end
+    end
+end