Mercurial > repos > public > sbplib

comparison +sbp/+implementations/d2_blocknorm_4.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.751e3 / 0.3456e4 0.661e3 / 0.3456e4 -0.515e3 / 0.3456e4 0.5e1 / 0.128e3; 0.661e3 / 0.3456e4 0.1405e4 / 0.1152e4 -0.3e1 / 0.128e3 0.29e2 / 0.3456e4; -0.515e3 / 0.3456e4 -0.3e1 / 0.128e3 0.989e3 / 0.1152e4 0.149e3 / 0.3456e4; 0.5e1 / 0.128e3 0.29e2 / 0.3456e4 0.149e3 / 0.3456e4 0.3407e4 / 0.3456e4;]; 8 H_U=[0.751e3 / 0.3456e4 0.661e3 / 0.3456e4 -0.515e3 / 0.3456e4 0.5e1 / 0.128e3; 0.661e3 / 0.3456e4 0.1405e4 / 0.1152e4 -0.3e1 / 0.128e3 0.29e2 / 0.3456e4; -0.515e3 / 0.3456e4 -0.3e1 / 0.128e3 0.989e3 / 0.1152e4 0.149e3 / 0.3456e4; 0.5e1 / 0.128e3 0.29e2 / 0.3456e4 0.149e3 / 0.3456e4 0.3407e4 / 0.3456e4;];
9 9
10 10
11 H=eye(m); 11 H=speye(m);
12 H(1:4,1:4)=H_U; 12 H(1:4,1:4)=H_U;
13 H(m-3:m,m-3:m)=flipud( fliplr(H_U(1:4,1:4) ) ); 13 H(m-3:m,m-3:m)=rot90( H_U(1:4,1:4) ,2 );
14 H=H*h; 14 H=H*h;
15 HI=inv(H); 15 HI=inv(H);
16 16
17 Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2)); 17 e=ones(m,1);
18 Q=spdiags([e -8*e 0*e 8*e -e], -2:2, m, m)/12;
19 % Q=(-1/12*diag(ones(m-2,1),2)+8/12*diag(ones(m-1,1),1)-8/12*diag(ones(m-1,1),-1)+1/12*diag(ones(m-2,1),-2));
18 20
19 Q_U = [-0.1e1 / 0.2e1 0.55e2 / 0.72e2 -0.47e2 / 0.144e3 0.1e1 / 0.16e2; -0.55e2 / 0.72e2 0 0.43e2 / 0.48e2 -0.19e2 / 0.144e3; 0.47e2 / 0.144e3 -0.43e2 / 0.48e2 0 0.47e2 / 0.72e2; -0.1e1 / 0.16e2 0.19e2 / 0.144e3 -0.47e2 / 0.72e2 0;]; 21 Q_U = [-0.1e1 / 0.2e1 0.55e2 / 0.72e2 -0.47e2 / 0.144e3 0.1e1 / 0.16e2; -0.55e2 / 0.72e2 0 0.43e2 / 0.48e2 -0.19e2 / 0.144e3; 0.47e2 / 0.144e3 -0.43e2 / 0.48e2 0 0.47e2 / 0.72e2; -0.1e1 / 0.16e2 0.19e2 / 0.144e3 -0.47e2 / 0.72e2 0;];
20 22
21 Q(1:4,1:4)=Q_U; 23 Q(1:4,1:4)=Q_U;
22 Q(m-3:m,m-3:m)=flipud( fliplr(-Q_U(1:4,1:4) ) ); 24 Q(m-3:m,m-3:m)=rot90( -Q_U(1:4,1:4) ,2 );
23 25
24 D1=HI*Q; 26 D1=H\Q;
25 27
26 M_U=[0.359e3 / 0.288e3 -0.443e3 / 0.288e3 0.97e2 / 0.288e3 -0.13e2 / 0.288e3; -0.51e2 / 0.32e2 0.325e3 / 0.96e2 -0.191e3 / 0.96e2 0.19e2 / 0.96e2; 0.43e2 / 0.96e2 -0.69e2 / 0.32e2 0.293e3 / 0.96e2 -0.137e3 / 0.96e2; -0.29e2 / 0.288e3 0.89e2 / 0.288e3 -0.427e3 / 0.288e3 0.727e3 / 0.288e3;]; 28 M_U=[0.359e3 / 0.288e3 -0.443e3 / 0.288e3 0.97e2 / 0.288e3 -0.13e2 / 0.288e3; -0.51e2 / 0.32e2 0.325e3 / 0.96e2 -0.191e3 / 0.96e2 0.19e2 / 0.96e2; 0.43e2 / 0.96e2 -0.69e2 / 0.32e2 0.293e3 / 0.96e2 -0.137e3 / 0.96e2; -0.29e2 / 0.288e3 0.89e2 / 0.288e3 -0.427e3 / 0.288e3 0.727e3 / 0.288e3;];
27 29
28 30
29 31
30 M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0)); 32 % M=-(-1/12*diag(ones(m-2,1),2)+16/12*diag(ones(m-1,1),1)+16/12*diag(ones(m-1,1),-1)-1/12*diag(ones(m-2,1),-2)-30/12*diag(ones(m,1),0));
33 M=-spdiags([-e 16*e -30*e 16*e -e], -2:2, m, m)/12;
31 34
32 M(1:4,1:4)=M_U; 35 M(1:4,1:4)=M_U;
33 36
34 M(m-3:m,m-3:m)=flipud( fliplr( M_U ) ); 37 M(m-3:m,m-3:m)=rot90( M_U ,2 );
35 M=M/h; 38 M=M/h;
36 39
37 DS_U=[0.25e2 / 0.12e2 -4 3 -0.4e1 / 0.3e1 0.1e1 / 0.4e1;]; 40 DS_U=[0.25e2 / 0.12e2 -4 3 -0.4e1 / 0.3e1 0.1e1 / 0.4e1;];
38 DS=zeros(m,m); 41 DS=sparse(m,m);
39 DS(1,1:5)=DS_U; 42 DS(1,1:5)=DS_U;
40 DS(m,m-4:m)=fliplr(DS_U); 43 DS(m,m-4:m)=fliplr(DS_U);
41 DS=DS/h; 44 DS=DS/h;
42 45
43 D2=HI*(-M+DS); 46 D2=H\(-M+DS);
44 47
45 d3=[-1 3 -3 1]; 48 % d3=[-1 3 -3 1];
46 t3=sum(abs(d3)); 49 % t3=sum(abs(d3));
47 DD_3(1:1,1:4)=[d3]; 50 % DD_3(1:1,1:4)=[d3];
48 DD_3(m:m,m-3:m)=[d3]; 51 % DD_3(m:m,m-3:m)=[d3];
49 52
50 % This works for wave eq. 53 % This works for wave eq.
51 % For studs interface in 1D no AD is needed. 54 % For studs interface in 1D no AD is needed.
52 ADD=1*h/(t3)*DD_3'*DD_3; 55 % ADD=1*h/(t3)*DD_3'*DD_3;
53 56
54 e_1 = zeros(m,1); 57 e_1 = sparse(m,1);
55 e_1(1)= 1; 58 e_1(1)= 1;
56 e_m = zeros(m,1); 59 e_m = sparse(m,1);
57 e_m(end)= 1; 60 e_m(end)= 1;
58 S_1 = -DS(1,:)'; 61 S_1 = -DS(1,:)';
59 S_m = DS(end,:)'; 62 S_m = DS(end,:)';
60 63
61 Q = H*D1-(-e_1*e_1' + e_m*e_m'); 64 Q = H*D1-(-(e_1*e_1') + (e_m*e_m'));
62 M = -(H*D2-(-e_1*S_1' + e_m*S_m')); 65 M = -(H*D2-(-e_1*S_1' + e_m*S_m'));
63 end 66 end