Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_blocknorm_6.m @ 267:f7ac3cd6eeaa operator_remake

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Sparsified all implementation files, removed all matlab warnings, fixed small bugs on minimum grid points.
author Martin Almquist <martin.almquist@it.uu.se>
date Fri, 09 Sep 2016 14:53:41 +0200
parents bfa130b7abf6
children
comparison
equal deleted inserted replaced
266:bfa130b7abf6 267:f7ac3cd6eeaa
6 end 6 end
7 7
8 H_U=[0.8489084265971e13 / 0.45952647390720e14 0.24636450459943e14 / 0.98469958694400e14 -0.2796787072531e13 / 0.12308744836800e14 0.2793599068823e13 / 0.14360202309600e14 -0.66344569931569e14 / 0.689289710860800e15 0.3784697867191e13 / 0.137857942172160e15; 0.24636450459943e14 / 0.98469958694400e14 0.27815394775103e14 / 0.19693991738880e14 -0.445601472229e12 / 0.861612138576e12 0.3896159037731e13 / 0.17232242771520e14 -0.866505556741e12 / 0.27571588434432e14 -0.25625418493681e14 / 0.689289710860800e15; -0.2796787072531e13 / 0.12308744836800e14 -0.445601472229e12 / 0.861612138576e12 0.31409405327129e14 / 0.17232242771520e14 -0.1595539040819e13 / 0.3446448554304e13 0.2651608170899e13 / 0.17232242771520e14 0.1434714163381e13 / 0.43080606928800e14; 0.2793599068823e13 / 0.14360202309600e14 0.3896159037731e13 / 0.17232242771520e14 -0.1595539040819e13 / 0.3446448554304e13 0.6984350202787e13 / 0.5744080923840e13 -0.62662743973e11 / 0.861612138576e12 -0.435331581619e12 / 0.12308744836800e14; -0.66344569931569e14 / 0.689289710860800e15 -0.866505556741e12 / 0.27571588434432e14 0.2651608170899e13 / 0.17232242771520e14 -0.62662743973e11 / 0.861612138576e12 0.20320736807807e14 / 0.19693991738880e14 0.1368363924007e13 / 0.98469958694400e14; 0.3784697867191e13 / 0.137857942172160e15 -0.25625418493681e14 / 0.689289710860800e15 0.1434714163381e13 / 0.43080606928800e14 -0.435331581619e12 / 0.12308744836800e14 0.1368363924007e13 / 0.98469958694400e14 0.27414523542149e14 / 0.27571588434432e14;]; 8 H_U=[0.8489084265971e13 / 0.45952647390720e14 0.24636450459943e14 / 0.98469958694400e14 -0.2796787072531e13 / 0.12308744836800e14 0.2793599068823e13 / 0.14360202309600e14 -0.66344569931569e14 / 0.689289710860800e15 0.3784697867191e13 / 0.137857942172160e15; 0.24636450459943e14 / 0.98469958694400e14 0.27815394775103e14 / 0.19693991738880e14 -0.445601472229e12 / 0.861612138576e12 0.3896159037731e13 / 0.17232242771520e14 -0.866505556741e12 / 0.27571588434432e14 -0.25625418493681e14 / 0.689289710860800e15; -0.2796787072531e13 / 0.12308744836800e14 -0.445601472229e12 / 0.861612138576e12 0.31409405327129e14 / 0.17232242771520e14 -0.1595539040819e13 / 0.3446448554304e13 0.2651608170899e13 / 0.17232242771520e14 0.1434714163381e13 / 0.43080606928800e14; 0.2793599068823e13 / 0.14360202309600e14 0.3896159037731e13 / 0.17232242771520e14 -0.1595539040819e13 / 0.3446448554304e13 0.6984350202787e13 / 0.5744080923840e13 -0.62662743973e11 / 0.861612138576e12 -0.435331581619e12 / 0.12308744836800e14; -0.66344569931569e14 / 0.689289710860800e15 -0.866505556741e12 / 0.27571588434432e14 0.2651608170899e13 / 0.17232242771520e14 -0.62662743973e11 / 0.861612138576e12 0.20320736807807e14 / 0.19693991738880e14 0.1368363924007e13 / 0.98469958694400e14; 0.3784697867191e13 / 0.137857942172160e15 -0.25625418493681e14 / 0.689289710860800e15 0.1434714163381e13 / 0.43080606928800e14 -0.435331581619e12 / 0.12308744836800e14 0.1368363924007e13 / 0.98469958694400e14 0.27414523542149e14 / 0.27571588434432e14;];
9 9
10 10
11 H=eye(m); 11 H=speye(m);
12 H(1:6,1:6)=H_U; 12 H(1:6,1:6)=H_U;
13 H(m-5:m,m-5:m)=flipud( fliplr(H_U(1:6,1:6) ) ); 13 H(m-5:m,m-5:m)=rot90( H_U(1:6,1:6) ,2 );
14 H=H*h; 14 H=H*h;
15 HI=inv(H); 15 HI=inv(H);
16 16
17 Q=(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)); 17 diags = -3:3;
18 stencil = [-1/60,3/20,-3/4,0,3/4,-3/20,1/60];
19 Q = stripeMatrix(stencil, diags, m);
20
21 % Q=(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));
18 22
19 Q_U = [-0.1e1 / 0.2e1 0.151864337282617e15 / 0.172322427715200e15 -0.251539972254817e15 / 0.344644855430400e15 0.61230525943549e14 / 0.114881618476800e15 -0.80987306509439e14 / 0.344644855430400e15 0.697178163343e12 / 0.13785794217216e14; -0.151864337282617e15 / 0.172322427715200e15 0 0.12350422095979e14 / 0.7658774565120e13 -0.78802251164141e14 / 0.68928971086080e14 0.4229407848431e13 / 0.7658774565120e13 -0.5372490790279e13 / 0.38293872825600e14; 0.251539972254817e15 / 0.344644855430400e15 -0.12350422095979e14 / 0.7658774565120e13 0 0.2217674201683e13 / 0.1723224277152e13 -0.13219134462287e14 / 0.22976323695360e14 0.19660399553981e14 / 0.114881618476800e15; -0.61230525943549e14 / 0.114881618476800e15 0.78802251164141e14 / 0.68928971086080e14 -0.2217674201683e13 / 0.1723224277152e13 0 0.62307836637379e14 / 0.68928971086080e14 -0.84068101764193e14 / 0.344644855430400e15; 0.80987306509439e14 / 0.344644855430400e15 -0.4229407848431e13 / 0.7658774565120e13 0.13219134462287e14 / 0.22976323695360e14 -0.62307836637379e14 / 0.68928971086080e14 0 0.44756810052211e14 / 0.57440809238400e14; -0.697178163343e12 / 0.13785794217216e14 0.5372490790279e13 / 0.38293872825600e14 -0.19660399553981e14 / 0.114881618476800e15 0.84068101764193e14 / 0.344644855430400e15 -0.44756810052211e14 / 0.57440809238400e14 0;]; 23 Q_U = [-0.1e1 / 0.2e1 0.151864337282617e15 / 0.172322427715200e15 -0.251539972254817e15 / 0.344644855430400e15 0.61230525943549e14 / 0.114881618476800e15 -0.80987306509439e14 / 0.344644855430400e15 0.697178163343e12 / 0.13785794217216e14; -0.151864337282617e15 / 0.172322427715200e15 0 0.12350422095979e14 / 0.7658774565120e13 -0.78802251164141e14 / 0.68928971086080e14 0.4229407848431e13 / 0.7658774565120e13 -0.5372490790279e13 / 0.38293872825600e14; 0.251539972254817e15 / 0.344644855430400e15 -0.12350422095979e14 / 0.7658774565120e13 0 0.2217674201683e13 / 0.1723224277152e13 -0.13219134462287e14 / 0.22976323695360e14 0.19660399553981e14 / 0.114881618476800e15; -0.61230525943549e14 / 0.114881618476800e15 0.78802251164141e14 / 0.68928971086080e14 -0.2217674201683e13 / 0.1723224277152e13 0 0.62307836637379e14 / 0.68928971086080e14 -0.84068101764193e14 / 0.344644855430400e15; 0.80987306509439e14 / 0.344644855430400e15 -0.4229407848431e13 / 0.7658774565120e13 0.13219134462287e14 / 0.22976323695360e14 -0.62307836637379e14 / 0.68928971086080e14 0 0.44756810052211e14 / 0.57440809238400e14; -0.697178163343e12 / 0.13785794217216e14 0.5372490790279e13 / 0.38293872825600e14 -0.19660399553981e14 / 0.114881618476800e15 0.84068101764193e14 / 0.344644855430400e15 -0.44756810052211e14 / 0.57440809238400e14 0;];
20 24
21 Q(1:6,1:6)=Q_U; 25 Q(1:6,1:6)=Q_U;
22 Q(m-5:m,m-5:m)=flipud( fliplr(-Q_U(1:6,1:6) ) ); 26 Q(m-5:m,m-5:m)=rot90( -Q_U(1:6,1:6) ,2 );
23 27
24 D1=HI*Q; 28 D1=H\Q;
25 29
26 M_U=[0.960901171090739e15 / 0.689289710860800e15 -0.502032138770899e15 / 0.229763236953600e15 0.493085196645929e15 / 0.344644855430400e15 -0.329491854944251e15 / 0.344644855430400e15 0.89541920186441e14 / 0.229763236953600e15 -0.50617198740721e14 / 0.689289710860800e15; -0.100483015499831e15 / 0.45952647390720e14 0.807564929223191e15 / 0.137857942172160e15 -0.415779274818991e15 / 0.68928971086080e14 0.80693719872887e14 / 0.22976323695360e14 -0.196663473955997e15 / 0.137857942172160e15 0.37943821632959e14 / 0.137857942172160e15; 0.99938177941669e14 / 0.68928971086080e14 -0.84419552767043e14 / 0.13785794217216e14 0.106922123424097e15 / 0.11488161847680e14 -0.223356054245897e15 / 0.34464485543040e14 0.157526160982357e15 / 0.68928971086080e14 -0.10062402380533e14 / 0.22976323695360e14; -0.68310884976863e14 / 0.68928971086080e14 0.17038649985979e14 / 0.4595264739072e13 -0.231397767539273e15 / 0.34464485543040e14 0.232669188399619e15 / 0.34464485543040e14 -0.1657930371065e13 / 0.510584971008e12 0.34774771016773e14 / 0.68928971086080e14; 0.18789143112277e14 / 0.45952647390720e14 -0.213895716727517e15 / 0.137857942172160e15 0.171024751153381e15 / 0.68928971086080e14 -0.8523669967037e13 / 0.2552924855040e13 0.485768751245399e15 / 0.137857942172160e15 -0.229158724354277e15 / 0.137857942172160e15; -0.51766014925489e14 / 0.689289710860800e15 0.202930494289627e15 / 0.689289710860800e15 -0.54332868549353e14 / 0.114881618476800e15 0.180479548146281e15 / 0.344644855430400e15 -0.1146942437956153e16 / 0.689289710860800e15 0.211001773091419e15 / 0.76587745651200e14;]; 30 M_U=[0.960901171090739e15 / 0.689289710860800e15 -0.502032138770899e15 / 0.229763236953600e15 0.493085196645929e15 / 0.344644855430400e15 -0.329491854944251e15 / 0.344644855430400e15 0.89541920186441e14 / 0.229763236953600e15 -0.50617198740721e14 / 0.689289710860800e15; -0.100483015499831e15 / 0.45952647390720e14 0.807564929223191e15 / 0.137857942172160e15 -0.415779274818991e15 / 0.68928971086080e14 0.80693719872887e14 / 0.22976323695360e14 -0.196663473955997e15 / 0.137857942172160e15 0.37943821632959e14 / 0.137857942172160e15; 0.99938177941669e14 / 0.68928971086080e14 -0.84419552767043e14 / 0.13785794217216e14 0.106922123424097e15 / 0.11488161847680e14 -0.223356054245897e15 / 0.34464485543040e14 0.157526160982357e15 / 0.68928971086080e14 -0.10062402380533e14 / 0.22976323695360e14; -0.68310884976863e14 / 0.68928971086080e14 0.17038649985979e14 / 0.4595264739072e13 -0.231397767539273e15 / 0.34464485543040e14 0.232669188399619e15 / 0.34464485543040e14 -0.1657930371065e13 / 0.510584971008e12 0.34774771016773e14 / 0.68928971086080e14; 0.18789143112277e14 / 0.45952647390720e14 -0.213895716727517e15 / 0.137857942172160e15 0.171024751153381e15 / 0.68928971086080e14 -0.8523669967037e13 / 0.2552924855040e13 0.485768751245399e15 / 0.137857942172160e15 -0.229158724354277e15 / 0.137857942172160e15; -0.51766014925489e14 / 0.689289710860800e15 0.202930494289627e15 / 0.689289710860800e15 -0.54332868549353e14 / 0.114881618476800e15 0.180479548146281e15 / 0.344644855430400e15 -0.1146942437956153e16 / 0.689289710860800e15 0.211001773091419e15 / 0.76587745651200e14;];
27 31
28 32
29 33
30 M=-(2*diag(ones(m-3,1),3)-27*diag(ones(m-2,1),2)+270*diag(ones(m-1,1),1)+270*diag(ones(m-1,1),-1)-27*diag(ones(m-2,1),-2)+2*diag(ones(m-3,1),-3)-490*diag(ones(m,1),0))/180; 34 % M=-(2*diag(ones(m-3,1),3)-27*diag(ones(m-2,1),2)+270*diag(ones(m-1,1),1)+270*diag(ones(m-1,1),-1)-27*diag(ones(m-2,1),-2)+2*diag(ones(m-3,1),-3)-490*diag(ones(m,1),0))/180;
35 diags = -3:3;
36 stencil = -1/180*[2,-27,270,-490,270,-27,2];
37 M = stripeMatrix(stencil, diags, m);
31 38
32 M(1:6,1:6)=M_U; 39 M(1:6,1:6)=M_U;
33 40
34 M(m-5:m,m-5:m)=flipud( fliplr( M_U ) ); 41 M(m-5:m,m-5:m)=rot90( M_U ,2 );
35 M=M/h; 42 M=M/h;
36 43
37 DS_U=[0.137e3 / 0.60e2 -5 5 -0.10e2 / 0.3e1 0.5e1 / 0.4e1 -0.1e1 / 0.5e1;]; 44 DS_U=[0.137e3 / 0.60e2 -5 5 -0.10e2 / 0.3e1 0.5e1 / 0.4e1 -0.1e1 / 0.5e1;];
38 DS=zeros(m,m); 45 DS=sparse(m,m);
39 DS(1,1:6)=DS_U; 46 DS(1,1:6)=DS_U;
40 DS(m,m-5:m)=fliplr(DS_U); 47 DS(m,m-5:m)=fliplr(DS_U);
41 DS=DS/h; 48 DS=DS/h;
42 49
43 D2=HI*(-M+DS); 50 D2=H\(-M+DS);
44 51
45 d5=[-1 5 -10 10 -5 1]; 52 % d5=[-1 5 -10 10 -5 1];
46 t5=sum(abs(d5)); 53 % t5=sum(abs(d5));
47 DD_5(1:1,1:6)=[d5]; 54 % DD_5(1:1,1:6)=[d5];
48 DD_5(m:m,m-5:m)=[d5]; 55 % DD_5(m:m,m-5:m)=[d5];
56 %
57 % % This works for wave eq.
58 % % For studs interface in 1D no AD is needed.
59 % ADD=7*h/(t5)*DD_5'*DD_5;
49 60
50 % This works for wave eq. 61 e_1 = sparse(m,1);
51 % For studs interface in 1D no AD is needed.
52 ADD=7*h/(t5)*DD_5'*DD_5;
53
54 e_1 = zeros(m,1);
55 e_1(1)= 1; 62 e_1(1)= 1;
56 e_m = zeros(m,1); 63 e_m = sparse(m,1);
57 e_m(end)= 1; 64 e_m(end)= 1;
58 S_1 = -DS(1,:)'; 65 S_1 = -DS(1,:)';
59 S_m = DS(end,:)'; 66 S_m = DS(end,:)';
60 67
61 Q = H*D1-(-e_1*e_1' + e_m*e_m'); 68 Q = H*D1-(-(e_1*e_1') + (e_m*e_m'));
62 M = -(H*D2-(-e_1*S_1' + e_m*S_m')); 69 M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
63 end 70 end