Mercurial > repos > public > sbplib

changeset 199:d18096820ed4 feature/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merged default into feature/grids.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 13 Jun 2016 16:50:16 +0200
parents 6fb354955c37 (current diff) 4f7930d2d2c4 (diff)
children ef41fde95ac4 38f203f00f3a
files
diffstat 8 files changed, 146 insertions(+), 3 deletions(-) [+]
line wrap: on
line diff
diff -r 6fb354955c37 -r d18096820ed4 +grid/bspline.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/+grid/bspline.m	Mon Jun 13 16:50:16 2016 +0200
@@ -0,0 +1,47 @@
+% Calculates a D dimensional p-order bspline at t with knots T and control points P.
+%  T = [t0 t1 t2 ... tm] is a 1 x (m+1) vector with non-decresing elements and t0 = 0 tm = 1.
+%  P = [P0 P1 P2 ... Pn] is a D x (n+1) matrix.
+
+% knots p+1 to m-p-1 are the internal knots
+
+% Implemented from: http://mathworld.wolfram.com/B-Spline.html
+function C = bspline(t,p,P,T)
+    m = length(T) - 1;
+    n = size(P,2) - 1;
+    D = size(P,1);
+
+    assert(p == m - n - 1);
+
+    C = zeros(D,length(t));
+
+    for i = 0:n
+        for k = 1:D
+            C(k,:) = C(k,:) + P(k,1+i)*B(i,p,t,T);
+        end
+    end
+
+    % Curve not defined for t = 1 ? Ugly fix:
+    I = find(t == 1);
+    C(:,I) = repmat(P(:,end),[1,length(I)]);
+end
+
+function o = B(i, j, t, T)
+    if j == 0
+        o = T(1+i) <= t & t < T(1+i+1);
+        return
+    end
+
+    if T(1+i+j)-T(1+i) ~= 0
+        a = (t-T(1+i))/(T(1+i+j)-T(1+i));
+    else
+        a = t*0;
+    end
+
+    if T(1+i+j+1)-T(1+i+1) ~= 0
+        b = (T(1+i+j+1)-t)/(T(1+i+j+1)-T(1+i+1));
+    else
+        b = t*0;
+    end
+
+    o = a.*B(i, j-1, t, T) + b.*B(i+1, j-1, t, T);
+end
\ No newline at end of file
diff -r 6fb354955c37 -r d18096820ed4 +util/calc_borrowing.m
--- a/+util/calc_borrowing.m	Mon Jun 13 16:49:11 2016 +0200
+++ b/+util/calc_borrowing.m	Mon Jun 13 16:50:16 2016 +0200
@@ -32,7 +32,7 @@
 
 
 %% 2nd order compatible
-[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher2_compatible(m,h);
+[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible2(m,h);
 S1 = S_1*S_1'  + S_m*S_m';
 S2 = S2_1*S2_1' + S2_m*S2_m';
 S3 = S3_1*S3_1' + S3_m*S3_m';
@@ -46,7 +46,7 @@
 
 
 %% 4th order compatible
-[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher4_compatible(m,h);
+[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible4(m,h);
 S1 = S_1*S_1'  + S_m*S_m';
 S2 = S2_1*S2_1' + S2_m*S2_m';
 S3 = S3_1*S3_1' + S3_m*S3_m';
@@ -59,7 +59,7 @@
 fprintf('\n')
 
 %% 6th order compatible
-[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher6_compatible(m,h);
+[H, HI, D1, D4, e_1, e_m, M4, Q, S2_1, S2_m, S3_1, S3_m, S_1, S_m] = sbp.higher_compatible6(m,h);
 S1 = S_1*S_1'  + S_m*S_m';
 S2 = S2_1*S2_1' + S2_m*S2_m';
 S3 = S3_1*S3_1' + S3_m*S3_m';
@@ -69,3 +69,27 @@
 fprintf('6th order compatible\n')
 fprintf('alpha_II:  %.10f\n',alpha_II)
 fprintf('alpha_III: %.10f\n',alpha_III)
+fprintf('\n')3
+
+
+
+
+
+% Ordinary
+
+for order = [2 4 6 8 10]
+    op = sbp.Ordinary(m,h, order);
+
+    S_1 = op.boundary.S_1;
+    S_m = op.boundary.S_m;
+
+    M = op.norms.M;
+
+    S1 = S_1*S_1'  + S_m*S_m';
+    alpha  = util.matrixborrow(M, h*S1);
+    fprintf('%dth order Ordinary\n', order)
+    fprintf('alpha:  %.10f\n', alpha)
+    fprintf('\n')
+end
+
+
diff -r 6fb354955c37 -r d18096820ed4 getVarname.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/getVarname.m	Mon Jun 13 16:50:16 2016 +0200
@@ -0,0 +1,7 @@
+function names = getVarname(varargin)
+    names = cell(size(varargin));
+
+    for i = 1:numel(varargin)
+        names{i} = inputname(i);
+    end
+end
\ No newline at end of file
diff -r 6fb354955c37 -r d18096820ed4 labelTask.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/labelTask.m	Mon Jun 13 16:50:16 2016 +0200
@@ -0,0 +1,12 @@
+function o = labelTask(in)
+
+    switch class(in)
+        case 'char'
+            fprintf(in);
+            o = tic();
+        case 'uint64'
+            o = toc(in);
+            fprintf(' - done %fs\n', o);
+        otherwise
+            error('Unknow input type: %s', class(in))
+    end
diff -r 6fb354955c37 -r d18096820ed4 prof.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/prof.m	Mon Jun 13 16:50:16 2016 +0200
@@ -0,0 +1,10 @@
+function prof(f)
+    profile on
+    try
+        f();
+        profile viewer
+    catch e
+        fprintf(2, '\n%s', getReport(e));
+        profile clear
+    end
+end
\ No newline at end of file
diff -r 6fb354955c37 -r d18096820ed4 saveFigurePosition.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/saveFigurePosition.m	Mon Jun 13 16:50:16 2016 +0200
@@ -0,0 +1,6 @@
+function saveFigurePosition()
+    defaultPosition = get(0,'defaultfigureposition');
+    f = gcf;
+    defaultPosition(1:2) = f.Position(1:2);
+    set(0,'defaultfigureposition',defaultPosition);
+end
\ No newline at end of file
diff -r 6fb354955c37 -r d18096820ed4 setFontSize.m
--- a/setFontSize.m	Mon Jun 13 16:49:11 2016 +0200
+++ b/setFontSize.m	Mon Jun 13 16:50:16 2016 +0200
@@ -1,3 +1,6 @@
+% setFontSize(fig, pts)
+% Sets all fontsizes within a figure to size pts
+% The default value for pts is 16.
 function setFontSize(fig, pts)
     default_arg('pts', 16);
     set(findall(fig,'-property','FontSize'),'FontSize',pts);
diff -r 6fb354955c37 -r d18096820ed4 time.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/time.m	Mon Jun 13 16:50:16 2016 +0200
@@ -0,0 +1,34 @@
+function t_out = time(f, n)
+    default_arg('n',1);
+
+    if n == 1
+        t = timeOnce(f);
+    else
+        t = timeAvg(f, n);
+    end
+
+    if nargout == 1
+        t_out = t;
+    else
+        fprintf('Elapsed time is %.6f seconds.\n', t)
+    end
+end
+
+function t = timeOnce(f)
+    s = tic();
+
+    f();
+
+    t = toc(s);
+end
+
+
+function t = timeAvg(f, n)
+    s = tic();
+
+    for i = 1:n
+        f();
+    end
+
+    t = toc(s)/n;
+end