Mercurial > repos > public > sbplib

changeset 319:fc91e3a2d44a feature/beams

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Cleaned d4_variable_6_min_boundary_points.
author Jonatan Werpers <jonatan@werpers.com>
date Mon, 26 Sep 2016 08:52:03 +0200
parents 99005a80b4c2
children 3955ccffb7c8
files +sbp/+implementations/d4_variable_6_min_boundary_points.m
diffstat 1 files changed, 30 insertions(+), 55 deletions(-) [+]
line wrap: on
line diff
--- a/+sbp/+implementations/d4_variable_6_min_boundary_points.m	Mon Sep 26 08:44:17 2016 +0200
+++ b/+sbp/+implementations/d4_variable_6_min_boundary_points.m	Mon Sep 26 08:52:03 2016 +0200
@@ -18,67 +18,44 @@
 
     % H?r med 6 RP ist?llet f?r 8 f?r D4 operatorn, dock samma randderivator
 
-    H = diag(ones(m,1),0);
-    H(1:6,1:6) = diag([13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200]);
-    H(m-5:m,m-5:m) = fliplr(flipud(diag([13649/43200,12013/8640, 2711/4320,5359/4320,7877/8640,43801/43200])));
-
-    e_1 = zeros(m,1);e_1(1)=1;
-    e_m = zeros(m,1);e_m(m)=1;
+    BP = 6;
+    if(m<2*BP)
+        error(['Operator requires at least ' num2str(2*BP) ' grid points']);
+    end
 
-    S_U = [-25/12, 4, -3, 4/3, -1/4]/h;
-    S_1 = zeros(1,m);
-    S_1(1:5) = S_U;
-    S_m = zeros(1,m);
-    S_m(m-4:m) = fliplr(-S_U);
-
-    S2_U = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
-    S2_1 = zeros(1,m);
-    S2_1(1:5) = S2_U;
-    S2_m = zeros(1,m);
-    S2_m(m-4:m) = fliplr(S2_U);
-
-    S3_U = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3;
-    S3_1 = zeros(1,m);
-    S3_1(1:5) = S3_U;
-    S3_m = zeros(1,m);
-    S3_m(m-4:m) = fliplr(-S3_U);
+    % Norm
+    Hv = ones(m,1);
+    Hv(1:6) = [13649/43200,12013/8640,2711/4320,5359/4320,7877/8640, 43801/43200];
+    Hv(m-5:m) = rot90(Hv(1:6),2);
+    Hv = h*Hv;
+    H = spdiag(Hv, 0);
+    HI = spdiag(1./Hv, 0);
 
 
-    x1=0.70127127127127;
+    % Boundary operators
+    e_l = sparse(m,1);
+    e_l(1) = 1;
+    e_r = rot90(e_l, 2);
 
-
-    D1=(1/60*diag(ones(m-3,1),3)-9/60*diag(ones(m-2,1),2)+45/60*diag(ones(m-1,1),1)-45/60*diag(ones(m-1,1),-1)+9/60*diag(ones(m-2,1),-2)-1/60*diag(ones(m-3,1),-3));
+    d1_l = sparse(m,1);
+    d1_l(1:5) = [-25/12, 4, -3, 4/3, -1/4]/h;
+    d1_r = -rot90(d1_l);
 
-    D1(1:6,1:9) = [
-        -21600/13649, 43200/13649*x1-7624/40947, -172800/13649*x1+   715489/81894, 259200/13649*x1-187917/13649, -172800/13649*   x1+735635/81894, 43200/13649*x1-89387/40947, 0, 0, 0;
-         -8640/12013*x1+7624/180195, 0, 86400/12013*x1-57139/12013,  -172800/12013*x1+745733/72078, 129600/12013*x1-91715/12013,  -34560/12013*x1+240569/120130, 0, 0, 0;
-          17280/2711*x1-715489/162660, -43200/2711*x1+57139/5422, 0,   86400/2711*x1-176839/8133, -86400/2711*x1+242111/10844,   25920/2711*x1-182261/27110, 0, 0, 0;
-         -25920/5359*x1+187917/53590, 86400/5359*x1-745733/64308,  -86400/5359*x1+176839/16077, 0, 43200/5359*x1-165041/32154,  -17280/5359*x1+710473/321540, 72/5359, 0, 0;
-          34560/7877*x1-147127/47262, -129600/7877*x1+91715/7877,   172800/7877*x1-242111/15754, -86400/7877*x1+165041/23631,   0, 8640/7877*x1, -1296/7877, 144/7877, 0;
-         -43200/43801*x1+89387/131403, 172800/43801*x1-240569/87602, -259200/43801*x1+182261/43801, 172800/43801*x1-710473/262806,  -43200/43801*x1, 0, 32400/43801, -6480/43801, 720/43801
-    ];
-    D1(m-5:m,m-8:m) = flipud( fliplr(-D1(1:6,1:9)));
-    D1 = D1/h;
+    d2_l = sparse(m,1);
+    d2_l(1:5) = [0.35e2/0.12e2 -0.26e2/0.3e1 0.19e2/0.2e1 -0.14e2/0.3e1 0.11e2/0.12e2;]/h^2;
+    d2_r = rot90(d2_l, 2);
 
+    d3_l = sparse(m,1);
+    d3_l(1:5) = [-0.5e1/0.2e1 9 -12 7 -0.3e1/0.2e1;]/h^3;
+    d3_r = -rot90(d3_l, 2);
 
-    H = h*H;
-    HI = inv(H);
-
-
-    D2 = [];
 
     % Fourth derivative, 1th order accurate at first 8 boundary points (still
     % yield 5th order convergence if stable: for example u_tt=-u_xxxx
 
-    m4 = 7/240;
-    m3 = -2/5;
-    m2 = 169/60;
-    m1 = -122/15;
-    m0 = 91/8;
-
-    M4=m4*(diag(ones(m-4,1),4)+diag(ones(m-4,1),-4))+m3*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3))+m2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2))+m1*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1))+m0*diag(ones(m,1),0);
-
-    %M4=(-1/6*(diag(ones(m-3,1),3)+diag(ones(m-3,1),-3) ) + 2*(diag(ones(m-2,1),2)+diag(ones(m-2,1),-2)) -13/2*(diag(ones(m-1,1),1)+diag(ones(m-1,1),-1)) + 28/3*diag(ones(m,1),0));
+    stencil = [7/240, -2/5, 169/60, -122/15, 91/8, -122/15, 169/60, -2/5, 7/240];
+    diags = -4:4;
+    M4 = stripeMatrix(stencil, diags, m);
 
     M4_U=[
         0.3504379e7/0.907200e6 -0.4613983e7/0.453600e6 0.4260437e7/0.453600e6 -0.418577e6/0.113400e6 0.524579e6/0.907200e6 0.535e3/0.18144e5;
@@ -90,10 +67,8 @@
     ];
 
     M4(1:6,1:6) = M4_U;
-
-    M4(m-5:m,m-5:m) = flipud( fliplr( M4_U ) );
-    M4 = M4/h^3;
+    M4(m-5:m,m-5:m) = rot90(M4_U, 2);
+    M4 = 1/h^3*M4;
 
-    D4 = HI*(M4-e_1*S3_1+e_m*S3_m  + S_1'*S2_1-S_m'*S2_m);
-
+    D4=HI*(M4 - e_l*d3_l'+e_r*d3_r' + d1_l*d2_l'-d1_r*d2_r');
 end